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ABSTRACT 



The application of analysis lattice filters to the 
problem of determining the input to a system from 
observations of the system's output (i.e., deconvolution) is 
discussed. Both linear and nonlinear systems are 
considered. Lattice filter modeling algorithms (Levinson 
and Schur) are presented. 

The theory of least-squares inverse filters is reviewed. 
This leads to a discussion of the lattice filter, which in 
turn leads to the Generalized Lattice Theory, The 
Generalized Lattice Theory is then used to develop a 
nonlinear lattice structure. Simulations show that the 
nonlinear lattice is an effective inverse filter for both 
linear and nonlinear systems. 
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I. INTRODUCTION 



T^he problem of estimating a signal based upon 
observations of a related signal is one of the most 
important operations in signal processing. The input and 
output signals of a linear system are related by the 
convolution operation 

y(t) = h(t) * x(t) =j h(t-T)x(T)d'r (1.1) 

— 00 

where h(t) is a causal, linear time-invariant (LTI), system 
impulse response. Since this thesis deals primarily with 
discrete digital signals, continuous time signals are 
sampled at uniform intervals, T, and are represented by 
discrete sequences x(nT) = x(n) for n = 0,1,2, ...,N. 

Discrete convolution for a LTI system is defined by 

n 

y(n) = h(n) * x(n) = h(n-m)x(m) (1*2) 

m=-oD 

and is shown in Figure 1.1. The basic deconvolution 
problem is to estimate the signal x(n), assuming that both 
y(n) and h(n) are known. Figure 1.2 depicts the inverse 
filtering process of recovering the input signal from the 
output signal by removing the system’s impulse response. 

Deconvolution has important applications in a variety 
of fields: Radar, communications, image processing, speech 
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Figure 1.1: Convolution 
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Figure 1.2 Deconvolution 



synthesis, seismology. For example, in image processing, 
deconvolution is used to recover the representation of the 
original object, x(m,n), from the representation of its 
image, y(m,n), by removing the blurring caused by the opti- 
cal system’s point-spread function, h(m,n). A common 
problem which arises in geophysics involves the deconvolu- 
tion of a seismic trace, y(n), into the approximately known 
impulsive waveform, h(n), and the desired reflection 
response, x(n), which reveals the structure of the layered 
Earth. In other applications, h(n) may represent the 
impulse response of a transmission channel, magnetic 
recording medium, or measurement device which broadens and 
smears (intersyrabol interference) the desired message x(n). 

There have been numerous approaches to the linear 
deconvolution problem, including least-squares filtering, 
linear inverse theory, linear programming, and homomorphic 
signal processing. The first part of this thesis deals with 
least-squares filtering techniques; the application of 
Kalman, waveshaping, and lattice filters to deconvolution is 
reviewed. Next, the lattice filter discussion is extended 
to the theory of the generalized lattice filter. Finally, 
nonlinear system theory is briefly reviewed, and the 
nonlinear lattice filter is developed and applied to the 
inverse filtering problem. Computer simulation results for 
both the linear and nonlinear lattice filters are presented. 
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II. LINEAR LEAST-SQUARES DECONVOLUTION TECHNIQUES 



A. INTRODUCTION 

This chapter begins with a discussion of the principles 
of least-squares filtering theory since deconvolution is 
essentially an inverse filtering process. Three specific 
types of least-squares filters are then introduced: Kalman, 
waveshaping, and lattice filters. The application of each 
of these three filter types to the problem of linear decon- 
volution is studied. Finally, the chapter concludes with 
the presentation of computer simulation results obtained for 
deconvolution experiments using the lattice filter. 

The goal of deconvolution is to recover the input 
signal, x(n), to a system based upon observed values of the 
system’s output signal y(n). An optimal processor must be 
determined to produce the best possible estimate of x(n) 
based upon present and past values of the output y<n). A 
traditional measure for defining the "best" signal processor 
is the minimum mean-squared error criterion. Using this 
criterion, the estimate x(n) is defined by a linear 
combination of the observed values y(n). Assuming 
causality, and that the observed signal is windowed to 
include only the M past samples, x(n) is written as 
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X ( n ) 



n > 0 



( 2 . 1 ) 



n 

y h ( n , i ) y ( i ) 
i=n-M 

where h(n,i) = 0 for n < i. The estimation error is given by 
e(n) = x(n) - x(n). To find the optimum signal processor, 
the h(n,i) coefficients which minimize the mean-square 
estimation error J, where 

2 2 

J = E[e (n)] = E[(x(n) - x(n)) ] , (2.2) 

must be determined. This is known as the Wiener filter 
formulation of the problem. [Ref. l:pp. 113,116] 

The least-squares theory of filtering began in the 
1940’s with the work of Norbert WIENER [Ref. 2:pp. 147-148]. 

WIENER developed a frequency domain procedure to design 
optimum filters, where optimality was defined by minimizing 
a mean-square error performance criterion. The Wiener 
filter is conventionally applied to linear time-invariant 
systems with stationary statistics when it is desired to 
separate one signal from another. In the early 1950’s, the 
Wiener filter was extended to include time varying and 
nonstationary statistics, but the calculations are 

cumbersome [Ref. 3:p. 1]. 

The mean-square estimation error J(n) is minimized by 
setting its partial derivatives, with respect to each of the 
filter coefficients h(n,i), equal to zero. Thus: 
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(3j(n) 

= 0 ( 2 . 3 ) 

C/h(n, i ) 



for i = n-M , n-M+ 1 , . . . , n . This yields a set of M+1 linear 
simultaneous equations, called the orthogonality equations 
where 



R (n,i) = E[e(n)y(i)] = 0 (2.4) 

ey 

for n-M i ^ n and n >_ 0 . R (n,i) is the 

ey 

crosscorrelation function between the error signal, e(n), 
and the data y(n). Substituting e(n) = x(n) - x(n) into the 
orthogonality equations gives the normal , or Wiener-Hopf 
equations : 

n 

E[x(n)y(i)] = V h ( n , k ) E [ y ( k ) y { i ) ] (2.5a) 

k = n-M 



or 

n 

R (n,i) = ^ h(n,k)R (k,i) . (2.5b) 

xy k=n-M yy 



Since the autocorrelation function R of the input signal 

yy 

and the crosscorrelation function R of the desired output 

xy 

signal with the input signal are known quantities, the M+1 
equations can be solved for the optimal filter weights 
h(n,i), i=n-M,...,n. If data vectors are defined so that 

T 

x(n) = [x(n-M), x ( n-M+ 1 ),..., x ( n ) ] (2.6a) 
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and 



T 

2 ;(n) = [y(n-M), y ( n-M+ 1 y ( n ) ] (2.6b) 

then the M+1 Wiener-Hopf equations can be written as 
T T 

E[x(n)y(n)]=hE[y(n)y;(n)] (2.7) 

where h = [h(n,n-M), h ( n-M+ 1 ),..., h ( n , n )] . Assuming that 
the signals x(n) and y(n) are stationary, then the filter 
coefficients are time-invariant and h(n,k) = h(n-k); and 

equation (2.7) can be written in matrix form as 



Ryy(O) Ryy(l) Ryy(2) ... Ryy(M) 




h(0) 




Rxy(O) 


Ryy(l) Ryy(O) Ryy( 1 ) ... 




h{ 1 ) 




Rxy( 1 ) 


Ryy(2) Ryy(l) Ryy(O) ... 

• • • • 




h( 2) 
• 


- 


Rxy(2) 

• 


« • • • 

• • • • 

Ryy(M) Ryy(O) 




• 

• 

h(M) 




• 

• 

Rxy(M) 



Now, the optimal coefficients can be solved by inverting 
the autocorrelation matrix, or by exploiting the matrix’s 
Toeplitz structure (all elements are the same on any given 
diagonal) to employ the more efficient Levinson algorithm. 
(Levinson’s algorithm will be discussed in the development 
of the analysis lattice filter.) 

The minimized value of the mean-square estimation error 
can now be computed, and is found to be 

2 n 

J (n) = E[x (n)] - ^ h(n,i )E[y( i )x(n) ] (2.9) 

min i=n-M 
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T T -1 

= R (0)-E[x(n)j; ( n ) ] E [ j^( n ) JL (n)] E [z( n ) x ( n ) ] 

XX 

[Ref. 4:p. 148]. These values correspond to the diagonal 

elements of the covariance matrix of the estimation error, 

T 

R = E[e(n)e (n)] (2.10) 

ee 

A 

where e(n) = x(n) - x(n) [Ref. l:p. 120]. Furthermore, the 
estimate of x(n) is given by 

A T -1 

x(n) = E[x(n)y (n)]E[y(n)y (n)] y(n) (2.11) 

= hy ( n) 

This estimate can be thought of as the projection of the 
desired signal x(n) onto the space spanned by the components 
of the observation vector y(n). The minimized estimation 
error vector is orthogonal, or normal, to the estimate 
x(n). [Ref. 4:p. 147] 

This completes the overview of least-squares filtering. 
The following sections present several linear deconvolution 
techniques which employ this criterion: Kalman filtering, 

spiking filters, and lattice filters. 

B. KALMAN FILTER 

In the early 1960 's, R.E. KALMAN introduced an optimal 
recursive filter based on state-space time-domain methods 
[Ref. 2;pp. 267-268]. The Kalman filter estimates the state 
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of a linear system, and is optimal in the sense that it min- 
imizes the mean-square error of the state estimate. The 
Kalman filter is useful when the system is defined by state 
space equations: The system signals are represented by 
random processes, and the data observations are 
contaminated by noise. The Kalman filter algorithm processes 
measurement data, and requires a priori state space models 
(known or assumed) of the system and measurement dynamics. 
Also, the statistics of the system input and measurement 
noises, as well as initial condition information, are 
required to produce the state estimate. This process is 
depicted in Figure 2.1. Here, the discrete Kalman filter 
equations will be presented, and then their application to 
deconvolution will be described. 

The state space representation of the discrete system 
and measurement models (see Figure 2.2) are written as [Ref. 



2 : pp . 195 


-200] : 






x(k) 


= F(k-l)x(k-l) + w(k-l) 


(2. 


, 12) 


^(k ) 


= H(k)x(k) + v(k) 


(2. 


, 13) 



where 



x(k) = 


(n 


X 


1 ) 


system state vector 


F(k) = 


(n 


X 


n) 


transition matrix 


w(k) = 


( n 


X 


1) 


system noise vector 


z( k ) = 


(m 


X 


1) 


measurement vector 
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Figure 2.1 



System, Measurement, and Kalman Filter 



x(t) = State Estimate 




^ I 

Figure 2.2 System and Measurement Model 
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H(k) 



(ra X n) observation matrix 
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v(k> = (m X 1) measurement noise vector. 

(Note that the time index k is used here to be consistent 

with the discrete Kalman filter literature.) 

The noise vectors w(k) and v(k) are assumed to have zero 

mean, white, Gaussian distributions with covariance matrices 

T T 

of Q(k) = E[w(k)w (k)] and R(k) = E[v(k)v (k)], 

respectively. Additionally, w and v are uncorrelated so 

T 

that E[w(k)v (j)] = 0 for all k and j . The state estimation 
error is defined by 

e(klk-l) = x(k) - x(k|k-l) (2.14) 

and the associated (n x n) error covariance matrix is 

T 

P(klk-l) = E[e(k|k-l)e (k|k-l)]. (2.15) 

An updated, a posteriori estimate of xd^l is obtained from 
the measurement z(k) and the a priori state estimate 
$ ( k I k- 1 ) by 

x(k) = x(k|k-l) + K( k ) [^( k ) -H ( k )^( k I k- 1 ) ] (2.16) 

where K(k) is the (n x m) Kalman gain matrix. The a 

posteriori error is given by 

e(k) = x(k) - x(k) (2.17) 

and the associated a posteriori error covariance matrix is 
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( 2 . 18 ) 



T 

P(k) = E[e(k)e (k) ] . 

The mean-square estimation error criterion is minimized 

when 

T T -1 

K(k)=P(k|k-l )H (k) [H(k)P(k|k-l )H (k)+ R(k)] (2.19) 

This optimal value of the Kalman gain matrix minimizes the 
individual terms along the main diagonal of P(k). 

Substituting the optimal gain matrix K(k) into the 

expression for P(k) results in 

P(k) = (i - K(k)H(k)) P(klk-l) , (2.20) 

where I_ is the identity matrix. In order to compute 
equations (2.16) and (2.19) recursively, the a priori 
estimates x(k+l|k) and P(k+l|k) must be determined at time 
k. The a priori estimates are given by 

x(k+l|k) = F(k)x(k) (2.21) 

T 

P(k+l|k) = E[e(k+l|k)e (k+l|k)] (2.22) 

T 

= E [ F ( k ) e ( k ) +w ( k ) ) ( F ( k ) e ( k ) +w( k ) ) ] 

T 

= F(k)P(k)F (k) + Q(k) . 

The Kalman filter algorithm is implemented by recursive- 
ly computing equations (2.16), (2.19), (2.20), (2.21), and 

(2.22). Figure 2.3 is a diagram of the Kalman filter. The 
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Figure 2.3 Discrete Kalman Filter 
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algorithm is initiated with the initial conditions 

x(0) = E[x(0) ] (2.23a) 

and , 

T 

P(0) = E[ (x(O)-x(O) )(x(0)-x(0) ) ]. (2.23b) 

In the event that a controlling input or a deterministic 
disturbance u(k) is applied to the system, the only change 
in the above algorithm is the state model and the a poster- 
iori estimate. They become [Ref. 5:p. 130] 

x(k) = F(k-l)x(k-l) + w(k-l) + G(k-l)u(k-l) (2.24) 

x(k) = F(k-1 )x(k-l )+G(k-l )u(k-l ) (2.25) 

+K(k) [z(k)-H(k)F(k-l )x(k-l ) ] . 

where G(k) is a (n x q) matrix and u(k) is a (q x 1) vector 
of input signals. 

The problem of estimating a desired signal based on 
noisy data observations pertains to such fields as 
communications, controls, and geophysics. However, the 
Kalman filter can also be applied to these deconvolution 
problems. As an example, the discrete Kalman inverse filter 
will now be applied to an exploration seismology problem. 

Typically, in the search for underground oil and gas 
deposits , a vibratory signal source generates a pulse of 
energy which is transmitted into the earth. In modern 
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seismology, the shape of this source wavelet can be 
carefully controlled; it is chosen so that it contains only 
those frequencies which are transmitted best by the earth. 
As the source wavelet propagates through the earth, it 
encounters many different layers with various acoustic 
impedances. At these layers, both partial reflection and 
refraction occur, creating numerous transmission paths. The 
received seismic signal at the surface is composed of many 
overlapping reflected wavelets. Therefore, the seismic 
trace can be represented as the convolution of the original 
source wavelet with an impulse train representing the 
various layers of the earth. Moreover, the seismic trace is 
contaminated by measurement noise and by the phenomena of 
ghost reflections and reverberations. [Ref. 6:pp. 14-15] 

The seismic trace can be described mathematically by 

z(t) = s(t,T)»r(t) + v(t) = / s(t,T)r(T)dT + v(t) (2.26) 

^ % 

where z(t) = measured seismic trace 

s(t,T) = finite duration, time varying wavelet 
r(t) = reflectivity function of earth’s structure 
v(t) = measurement noise. 

The seismologist must extract the structure of the earth, 
r(t), by analyzing the noisy seismic data z(t). This 
process of removing the wavelet shape from the trace and 
leaving behind the impulse train representing the reflected 
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wavelet’s strength and arrival time is that of 
deconvolution. CRUMP [Ref. 3:pp. 6-7] shows that if the 
seismic signal is sampled at discrete, uniformly spaced 
intervals, and if s(t,T) and r(t) are assumed to be causal, 
then 2 (t) is represented by 
J 

z(k) = Y [ s ( k , k-i+ 1 ) r ( k-i + 1 ) ] + v(k) (2.27) 

i = l 

where the sample number k = 1,2,3,..., and 

L = length of the wavelet given in number of samples 
J = k for k<L 
= L for k^L . 

When M traces of K samples in length are available for 
processing then z.(k) becomes a (M x 1) vector where the j-th 
component is given by 
J 

z (k) = Y. ( k) )r (k-i + 1 ) ] + v (k) (2.28) 

j i=l ji j 

for j = 1,2,...,M and k = 1,2,...,K. This assumes that the 
reflectivity function is the same for each trace, while the 
shape of the exciting wavelet may vary from trace to trace. 
The time-varying wavelet values are contained in the (M x L) 
matrix H(k) where the j-th row contains the L samples of the 
wavelet which generates the j-th trace: 

H (k) = s (k, k-i+1 ) . (2.29) 

ji j 

In vector form, the equations become 
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(2.30) 



z(k) = H(k)x(k) + v(k) 



where the state, measurement, and noise vectors are given by 

T 

x(k) = [r(k) ,r(k-l ) , . . . ,r(k-L+l ) ] (2.31a) 

T 

z(k) = [z (k),z (k),...,z (k)] (2.31b) 

1 2 M 

T 

v(k) = [v (k),v (k),...,v (k)] . (2.31c) 

12 M . 

This is the Kalman filter measurement model. Now the state 
model must be determined. 

The state model is arrived at by assuming a general 
relationship for the reflection coefficients of x<^). The 
assumed relationship is the autoregressive equation 
L 

r(k) = V [b (k-l)r(k-i)] + w(k-l) . (2.32) 

i = l i 

Comparing equation (2.32) with the state vector x ( k ) yields 
the state model 



x(k) = F(k,k-1 )x(k-l ) + w(k-l) (2.33) 

where the (M x L) transition matrix and the (M x 1) system 
noise vector are given by 



F(k,k-1) 



b (k-1 ) b (k-1 ) ... 


: b (k-1) 


1 2 


: L 


1 


! 0 




1 

1 



(2.34a) 
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w(k-l) 



T 



(2.34b) 



= [w(k-l ) , 0 , 0 , . . . , 0] , 

i is the identity matrix, and 0 is the null vector. 
Equations (2.30) and (2.32) provide the measurement and 
system models for implementing the deconvolution via the 
Kalman filter. CRUMP discusses methods by which to obtain 
numerical values for the reflection coefficient vector b(k) 
and the time-varying wavelet sample matrix H(k) [Ref. 3 : pp 
8-11]. Once these matrices are determined, the recursive 
Kalman filter removes the effects of s(t,T) from r(t) and 
generates the state estimate x(k) which provides L samples 
of the desired reflectivity function at each time k. 

It is interesting to note that both the Kalman and 
Wiener filters are minimum mean-square error estimators, 
both require the same a priori knowledge of the process to 
be estimated, and that both yield identical estimates. 
However, the Kalman filter does have distinct advantages 
over the Wiener filter. First, due to the matrix form of 
the state space equations, the Kalman filter has 
multichannel capability and is equivalent to a bank of 
optimal estimators. Moreover, the Kalman filter is ideally 
suited to computer implementation due to its discrete and 
recursive characteristics. [Ref.2:pp. 268-269] 

The discrete least-squares approach which follows is a 
viable alternative to the Kalman filter algorithm. 
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particularly when state space model equations are not 
available or applicable. 



C. LEAST-SQUARES INVERSE FILTER 

When a linear system, H(z), is excited by the an input 

signal x(n), the output y(n) is defined by the convolution 

relationship y(n) = x(n) * h(n). Deconvolution involves 

finding the inverse filter G(z) such that H(z)G(z) = 1. 

(Note that the discrete time domain is related to the fre- 

jwT 



quency domain through the equation z = e ' , where T is 
the discrete sampling interval.) This condition transforms 
to h(n) * g(n) = d(n) in the discrete time domain; h(n) and 
g(n) are the impulse responses of the filters H(z) and G(z), 
respectively, and d(n) is the unit impulse function. If 
this condition is met, then the original input signal x(n) 
is recovered at the output of the inverse filter. 

The transfer function of the causal system is defined by 
the infinite series obtained by taking the one-sided z- 
transform of the system’s impulse response: 

oo 



-1 



H(z) = Y(z)/X(z) = Y. 

i = 0 



(2.35) 



H(z) can be determined by inserting a known sequence x(n) 
into the system, measuring the output sequence y(n), and 
then manipulating their z-transforms . The inverse filter 
G(z) is then computed by carrying out the polynomial 
division 
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(2.36) 



-1 -2 

G(z) = 1/H(z) = l/[h(0)+h(l)z +h(2)z +...] 

-1 -2 -M 

= g(0) + g(l)z + g(2)z + ... + g(M)z + ... 

and truncating to M+1 terms if necessary. If the exact 
inverse filter is approximated by truncating G(z) to order 
M, then its impulse reponse is given by the sequence g(n), 
n=0 , 1 , 2 , . . . ,M . Furthermore, if the original system’s 

impulse response is represented by h(n), n=0 , 1 , 2 , . . . N , then 

M 

d(n) = g(n) * h(n) = V g(m)h(n-M) (2.37) 

m=0 

for 0 ^ n ^ N+M. This approximation to the impulse function 

improves as the order M of the inverse filter is increased. 

Now the stability of the inverse filter will be 

addressed. If H(z) has all its zeros inside the unit circle 

in the complex z-plane, it is referred to as a minimum-delay 

polynomial; the corresponding sequence h(n) is called a 

minimum phase-lag sequence. This is a sufficient condition 

to guarantee that H(z) has a stable inverse, because the 

zeros of H(z) become the poles of G(z), and G(z) is a stable 

filter if all its poles lie within the unit circle. 

Maximum- and mixed-delay signals are obtained by 

t 

transforming the zeros of H(z) from z to 1/z where the 

i i 

superscript represents the complex conjugate operation. 

A maximum-delay sequence has all of its zeros outside the 
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unit circle, while the mixed-delay sequence has zeros inside 
and outside the unit circle. 

If H(z) has N zeros, then transforming these zeros 

N 

results in at most 2 distinct sequences. Of note, each of 

these minimum, maximum, and mixed-delay signals have the 

2 * 

same magnitude spectrum: |H(w)| = H(w)H (w) [Ref. l:p. 98]. 

However, they do have distinct phase spectra [Ref. 4:p. 

175]. The maximum-delay polynomial can be written as 

R * t * -N 

H (z) = h (N) + h (N-l)z +...+ h (0)z (2.38) 

R 

The so called reverse polynomial H (z) is a conjugated, 

reflected, and shifted version of H(z). The corresponding 

R * * t 

maximum phase-lag sequence is h = {h (N),h (N-1 ),..., h (0)} 
[Ref. 6:p. 72] While the minimum-delay filter has a causal, 

stable inverse consisting only of a memory function, the 
maximum-delay filter has an inverse which consists only of a 
stable, noncausal, anticipation function. The stable 
inverse of a mixed-delay function consists of both memory 
and anticipation functions. Filters with nonvanishing 

anticipation components are noncausal; they cannot work in 
real time since the future values of the filter input are 
not available for processing. This problem can be 
circumvented if the entire signal is first recorded prior to 
analysis; then the required future input data is available. 
[Ref. 4:p. 87] 
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The energy distribution in minimum, mixed, and maximum 
delay signals will now be examined. Since each of these 
signals have an identical magnitude spectrum, they also 
have the same total energy. However, although the total 
energy is the same, the rate at which the energy builds up 
differs for the various sequences. Parseval’s theorem 
states that the total energy in a signal is given by 




IT 

|H(w) 



2 

dw/(2n) 



M 2 

Y, lh(m) I 
m=0 



(2.39) 



If the partial energy is defined as 
n 2 

P(n)=yih(m)|, (2.40) 

m=0 



then it can be shown that the energy builds up quickest in 
the minimum-delay sequence, and that it builds up the 
slowest in the maximum-delay sequence [Ref. 6:pp. 75-76]. In 
other words, the minimum-delay signal makes its impact as 
soon as possible since its energy is concentrated at the 
front of the sequence. The maximum-delay signal makes its 
major impact at a later time since its energy is 
concentrated at the end of the sequence. The energy curves 
associated with all the possible mixed-delay signals lie 
between these two extremes. Finally, it can be shown that 
the convolution of two minimum-delay sequences with one 
another results in a minimum-delay sequence. The convolution 
of maximum-delay signals results in a maximum-delay 
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sequence. Convolution involving any other combination of 
sequences results in a mixed-delay sequence. Of note, the 



resulting maximum and minimum-delay sequences are the 
reverse of each other. [Ref. 6:pp. 73-74] 

The method described above of finding an approximate, 
finite length inverse filter consisted of simply truncating 
the exact inverse filter found by polynomial division. It 
was seen that the inverse filter G(z) attempted to transform 
the impulse response of H(z) into a unit impulse located at 
the origin. This can be thought of as an attempt by G(z) to 
undo the blurring effect of H(z) (i.e., H(z) "blurs" the 
impulse d(n) into the impulse response h(n)). If the input 

to H(z) is designated x(n), and .if the output of G(z) is 
x(n), then the error of the approximated inverse filter is 

e(n) = x(n) - x(n) for 0 £ n ^ N+M. (2.41) 

The error energy is defined by 



For the polynomial division / truncation method, J decreases 
as the order M of G(z) increases [Ref. 6:p. 136]. Seeking to 
minimize the error energy J with respect to the inverse 
filter coefficients leads to another approach for finding an 
approximate inverse filter. 



J 




N+M 2 
y e (n) 



(2.42) 
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Filters which minimize the mean-square error J are 
called least error energy or least-squares inverses. In 
general, the desired output sequence x(n) of the inverse 
filter G(z) could be of any shape. G(z) is then called a 
waveshaping filter. When applied to deconvolution, the 
desired output of the inverse filter is a unit impulse 
d(n-i), where i = 0 , 1 , 2 , . . . , N+M defines the lag of the digi- 
tal inverse filter. In this case, G(z) is called an i-th 
delay spiking filter since it tries to condense the system 
impulse response h(n) into a spike with i delays. Since i = 
0 , 1 , 2 , . . . , N+M , there are N+M+1 possible spiking filters. As 
will be discussed, there are preferred values of the lag i 
for minimizing J, depending on the phase characteristics of 
h(n). The spiking filter problem is shown in Figure 2.4. 

It is convenient to restate the deconvolution problem in 
a matrix form based on the Yule-Walker, or autocorrelation, 
method [Ref. iJpp. 243-245]. This is accomplished by 
defining the (M+N+1 x M+1) system output matrix Y, the 
(M+1 X 1) impulse response vector g_, and the (N+M+1 x 1) 
input vector x as follows: 
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Figure 2.4 Spiking Filter 
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d(n) = unit impulse function 

h(n) = impulse response of system 

g(n) = impulse response of inverse filter 



< M ZEROS > 

y(0) 0 0 ... 0 

yU) y(0) 0 ... 0 

y(2) y(l) y(0) ... 0 



Y 



: : : ... o 

y(N) y(N-l) y(N-2) ... y(N-M) 

0 y(N) y(N-l) 

0 0 y(N) 

. . 0 



0 0 0 



y(N) 



(2.43a) 



T 

^ = [g(0), g(l), ... ,g(M)] , 



(2.43b) 



T 

X = [x(0) , x( 1 ) , ... ,x(N+M) ] 



(2.43c) 



The above notation is for the general case of the 

waveshaping filter. For the special case of the spiking 
filter, . the y(n) components of the system output matrix Y 

are replaced by the impulse response h(n). Also, the 
desired signal x reduces to an impulse d(n-i). Now 

equations (2.37), (2.41), and (2.42) can be rewritten as 

X = Yg: , (2.44a) 

e = X - X , and (2.44b) 

T 

J = e e . (2.44c) 



As discussed in the introduction to this chapter, the 
least-squares criterion is satisfied by minimizing J with 
respect to the inverse filter coefficients G. This results 
in the normal equations 
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(2.45) 



T T 

Y = Y X 



T 

which can be solved for g_. By recognizing that Y Y is 

equivalent to the (M+1 x M+1 ) sampled autocorrelation matrix 
T 

R = E[yx ] where the M+1 length data vector is = [ y(0), 

T T 

y ( 1 y ( N ), 0 , 0 0 ] , and that Y x is a length (M+1) 

cross-correlation column vector r, it can be seen that 



T -1 T -1 

g.= (YY) Yx = R r 


(2.46) 



-1 

where R can be evaluated efficiently by Levinson’s 



algorithm. The actual filter 


output is then 


-1 T 

^ = Yg. = (YR Y )x = Px 


(2.47) 


T -1 T 





where P = Y(Y Y) Y is a square (N+M+1) dimensioned matrix 
called the projection or performance matrix. The filter’s 
performance improves as P approaches the identity matrix. 
Now, the estimation error and cost function are written as 



e= x-^= x(I - P) 


(2.48) 


T T 

J = ee = X (I.~ P)x 


(2.49) 


Since in the deconvolution 


problem x is the spike with i 


delays, the inverse filter 


output X is actually the i-th 

-1 T 



column of the P matrix. Also, since = R Y , the 

coefficients of the i-th spiking filter are contained in the 

-1 T 

i-th column of the matrix R Y . Moreover, the energy 
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error of the i-th spiking filter reduces to J = 1 - P(i,i). 
Therefore, in order to realize the best inverse function 
(i.e., minimize J), select the delay i for which P(i,i) is 
largest. For a chosen lag i, the filter’s performance also 
improves as the order M of the inverse filter is increased. 

Now let J(i) represent the estimation error for the i-th 
spiking filter. The grand sum of squared errors is defined 
as V = J(0) + J(l) +...+ J(M+N). It can be shown that V = 
(M+N+1) - (M+1) = N, where N is the order of the system H(z) 
[Ref. 4:p. 198]. Therefore, V is independent of order M. 

For sufficiently long spiking filters, the optimal value 
of the delay i depends on the phase characteristics of the 
signal h(n) and the choice of the lag is governed by the 
following rules. If h(n) is a minimum-delay input, the 
spiking filter should have zero delay, i = 0. This says 
that for a signal with its energy concentrated towards the 
front of the sequence, it is easiest to condense it to a 
unit impulse at the origin. If the signal is a maximum- 
delay input, the maximum-delay spike i = N+M+1 should be 
selected. If h(n) is a mixed-delay signal, then the i 
corresponding to the largest P(i,i) should be chosen. [Ref 
6 . : p . 152] 

Under certain conditions, the estimation error will go 
to zero as the length of the inverse filter tends to 
infinity. First, as previously discussed, if h(n) is a 
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minimum-delay sequence then an exact inverse can be found 
through polynomial division. The zero-delay spiking filter 
approaches this exact inverse filter G(z) as the number of 
terms M+1 increases. Therefore, the estimation error goes 
to zero as M goes to infinity. The second condition is, if 
h(n) is not a minimum-delay sequence, J will approach zero 
as 1/M if the lag i of the spiking filter is chosen to be 
sufficiently large. [Ref.4:pp. 200-201] 

Up until now, the effects of measurement noise and 
imperfect knowledge of the distorting function H(z) have 
been neglected. If noise is introduced, then the output 
y(n) of the system H(z) driven by input signal x(n) becomes 

y(n) = h(n)*x(n) + v(n) (2.50) 

where v(n) is assumed to be zero-mean white noise with 
variance Q. If this signal is passed through the previously 
determined inverse filter G(z), the filter output becomes 

x(n) = g(n)*y(n) = g ( n ) *h ( n ) *x ( n ) + g(n)*v(n) (2.51) 

= d(n)*x(n) + u(n) ~ x(n) + u(n) 

where u(n) = g(n)*v(n) is the filtered noise signal and d(n) 
is the unit impulse function. The variance of u(n) is: 
2 M2 T 

E[u (n)] = Q 5] g (n) = Qg: g. (2.52) 

n = 0 

[Ref. l:p. 248]. This variance may be larger than the 
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original noise variance Q. For example, if h(n) is a low 
frequencey signal, then g(n) must be very spiky in order to 
compress h(n) into a high frequency content spike. As a 
result, g(n) will likely have values greater than one. 
Therefore, the variance of u(n) will be larger than Q. This 
further degrades the estimate x(n). 

To compensate for this filtered noise, the minimization 
criterion is modified so that 

N+M 2 T 

J = (d(n) - g(n)*h(n)) + AQg. g_ (2.53) 

n = 0 

where A is a positive parameter. The first term of the cost 
function tries to produce a good inverse filter whereas the 
second term tries to reduce the output noise. If A is large, 
noise reduction is emphasized at the expense of obtaining a 
good inverse function. A small A emphasizes finding a good 
inverse function g(n), and there is little output noise 
reduction . 

Using this new cost function, the resulting normal 
equations are given by 
T T 

(Y Y + AQDm. = Y X . (2.54) 

Comparing this with equation (2.45) reveals that the main 
diagonal elements of the sampled autocorrelation matrix R 
have been modified by an additive term: R(0) becomes R(0) + 

Aq. If the Backus-Gilbert "prewhitening" parameter epsilon 
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is introduced, where 



€ = AQ/R(0) , (2.55) 

then the diagonal elements are written as (1 +£ )R(0). Even 

very small values for the BG parameter have a beneficial 

T 

effect of stabilizing the inverse of the matrix (Y Y + AQX) • 
[Ref. l:p. 246] 

D. LATTICE FILTER 

1 . Introduction 

The lattice filter is another optimal least-squares 
predictor which can be applied to the linear deconvolution 
problem. The lattice, or ladder filter derives its name 
from the cascade form of its signal flow graph. Basically, 
the lattice filter provides an alternative to the 
transversal filter for modeling a signal. It can be viewed 
as a Gram-Schmidt orthogonalization of the incoming data. 
This will be elaborated upon in subsequent paragraphs. To 
date, much of the work done with lattice filters has been in 
the areas of speech analysis/synthesis, seismology, and in 
high-resolution spectral estimation [Ref. 7:p. 841]. Prior 
to developing the lattice filter equations, key mathematical 
concepts which apply to this discussion will be reviewed. 

2 . Mathematical Background 

A central concept behind the development of the 
lattice filter is that of orthogonality. Two vectors are 
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orthogonal if their inner product is zero. The geometric 
interpretation of this is that the vectors are at right 
angles to one another. As an example, if the random varia- 
bles u and V are the basis vectors of a two-dimensional 

space, they are orthogonal if and only if their inner 

product <u,v> = E[uv] = 0. In the case where the linear 
space is spanned by the random variables of a data vector, 
two vectors are orthogonal if the corresponding random 

variables are uncorrelated and if one or both have zero 
mean [Ref. 8:p. 92]. 

A related theorem is the orthogonal decomposition 

theorem, which states that any random variable may be 

decomposed uniquely with respect to a subspace S into two 

mutually orthogonal parts, one part which is parallel to S 

(i.e., lies in S) and the other part orthogonal to S. That 

IS, y can be written y = y + e where e is orthogonal to y 

and to the basis vectors which span the subspace S, and 

where y is defined by a linear combination of the basis 

vectors. If S is spanned by the basis vectors 

M 

{u ( 1 ), u ( 2 ),..., u ( M )} , then y = V a(i)u(i) where it can be 

i=l 

shown that the coefficients are given by the equation 

2 -1 

a(i) = E[yu(i)]E[u (i)j . [Ref. l:p. 13] 

The orthogonal projection theorem adds to this by 
stating; the orthogonal projection y of a vector y onto a 
linear subspace S is that vector in S which lies closest to 
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^ with respect to the distance induced by the inner product 
of the vector space [Ref, l:*p, 15], Restated, this simply 
means that y is the best estimate of y that can be made by a 
linear combination of the basis vectors of S in a minimum 
mean-squared error sense. Figure 2,5 shows the projection 
of y into the subspace S, 

Another important concept in developing the lattice 
filter is that of Gram-Schmidt orthogonalizat ion , The Gram- 
Schmidt procedure is a recursive orthogonalization process 
which generates a set of mutually orthogonal basis vectors 
{u( 1 ) ,u( 2 ) , , , , ,u(M) ) from a given set of basis vectors 
{y( 1 ) j y( 2 ),,,,, y(M) } , The procedure is initialized by 
letting u(l) = y(l). Next, y(2) is decomposed with respect 
to u(l). That part of y(2) which is orthogonal to u(l) 
becomes u(2). Next, y(3) is decomposed with respect to the 
subspace spanned by {u(l),u(2)). That part of y(3) which is 
orthogonal to this subspace becomes u(3). In general, the 
new set of orthogonal basis vectors is defined by 



u(n) r y(n) 



n- 1 



b( n, i )u( i ) 



(2.56) 



2 -1 

for 2 j< n ^ M and where b(n,i) = E[y(n)u(i)]E[u (i)] 

Using the orthogonal decomposition theorem, this is 
equivalent to y(n) = y(n) + u(n) where y(n) is the best 
estimate of the n-th component of the data vector based on 
the previous n-1 components and u(n) represents the 
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Figure 2.5 



Projection of y onto Space S 



41 



estimation error. If b(n,n) 



1 is introduced, then 



n 

y(n) = y b(n,i)u(i) 

i^l 



This can be written in matrix 



for 1 <. n ^ M. 
form as 



( 2 . 57a) 



Y_ = Bu where 

T 

z = [y( 1 ) . y( 2) , . . . , y(M) ] 

T 

u = [u( 1 ) , u( 2 ) , . . . , u(M) ] 



B 



1 0 0 

b( 2 , 1 ) 1 0 

b(3, 1) b(3,2) 1 

• • 

• • 

b(M,l) b(M,2) 



(2.57b) 



0 

0 

0 



1 



This is a convenient notation for representing the 

transformation from a set of correlated basis vectors y to 

an uncorrelated set of basis vectors u. The bases y and u 

span the same M-dimensional subspace S, but there are no 

redundant correlations between the basis vectors of set u. 

Since the basis of u is uncorrelated, its components u(i) 

(i = 1,2,...,M) are referred to as innovations because each 

additional component contributes completely new information 

to the estimate of y. [Ref. l:pp. 16-18] 

Finally, the transformation y = Bu corresponds to a 

LU (lower upper ) -Cholesky factorization of the correlation 

T 

matrix of y. By definition R = E[yy ]. Substituting for 

yy 
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Y. results in 



T T T 

R = BE[uu ]B = BR B (2.58) 

yy uu 

T 

where B is lower triangular, B upper triangular, and R is 

uu 

a diagonal matrix since the basis vectors u(i) are 

uncorrelated with one another. 

3 . Derivation of Lattice Filter Equations 

Now the lattice filter order update equations will 

be developed. A random signal y(n) can be modeled as the 

output of a causal, stable, linear filter H(z) which is 

driven by a stationary, uncorrelated, white noise sequence 

{e( 1 ), e( 2 ),..., e(P) ) [Ref. l:p. 30]. A P-th order 

autoregressive model is defined by H(z) = 1/A (z) where 

P 

-1 -2 -P 

A (z) = l+a(l)z +a(2)z +...+a(P)z . (2.59) 

P 

The filter A (z) has several names: prediction error filter, 

P 

inverse filter, or analysis filter. The signal y(n) can be 
written as 



y(n ) = e(n) - 
If the predictor y(n) 



P 




is 



a( i)y(n-i ) 
introduced , 



this becomes 



(2.60) 



y(n) - y(n) = e(n) (2.61) 

A P 

where y(n) = - ) a(i)y(n-i). Now, y(n) is the estimate of 

i = 1 

y(n) at time (n-1) based on the previous P samples, (y(n-P), 
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y ( n-P+ 1 ) , . . . , y(n-l)} and e(n) is equivalent to the estima- 

tion error. The estimate y(n) can be thought of as the 
projection of the ( P+ 1 ) -dimensional y(n) vector onto the P- 

dimensional subspace spanned by the components of the data 

T 

vector Y = [ y { n- 1 ), y ( n-2 y ( n-P ) ] . In general, the 

past values of y(n) are correlated with one another. There- 
fore, the random variable components of Y do not generally 
form a set of mutually orthogonal basis vectors. The Gram- 
Schmidt orthogonalization procedure removes these 
correlations and transforms Y into a set of mutually 
orthogonal basis vectors which span the same subspace. 
Additionally, the predictor coefficients a(i) are replaced 

by the reflection coefficients K . This results in the 

i 



lattice filter. 

In order to obtain the optimal estimator, the 

predictor coefficients which minimize the mean-squared 

2 

prediction error E[e (n)] must be determined. This problem 
was discussed in the introduction to this chapter. The 
resulting matrix form of the normal equations is 



R(0) R(l) R(2) . . . R(P) 




1 




Ep 


R( 1) R( 0 ) R( 1 ) ... 




a ( 1 ) 




0 


R(2) R(l) R(0) ... 




a ( 2 ) 


— 


0 


• • • • 

# • ♦ # 




• 

» 




• 

# 


R(P) R(0) 




a( P) 




0 



where R(k) = E [ y ( n+k ) y ( n ) ] and E is the minimized mean- 

P 

squared prediction error for the P-th order filter. The 
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sample autocorrelation matrix components are calculated from 
the length N data sequence by 

N-l-k 

R(k) = (1/N) y y(n+k)y{n) (2.63) 

n = 0 



for 0 ^ k <_ P and where P ^ N-1. [Ref. l:p. 150] 

Now there are P+1 equations and P+1 unknown model 

parameters { a ( 1 ) , a ( 2 ) , . . . , a ( P ) ; E }. The P+1 equations can be 

P 

solved by inverting the autocorrelation matrix. This 
3 2 

requires 0(P ) operations and 0(P ) storage locations. The 

P+1 equations can be solved more efficiently by taking 

advantage of the matrix’s Toeplitz structure by using 

Levinson’s algorithm. Levinson’s algorithm reduces the 

2 

required number of operations and storage locations to 0(P ) 
and 0(P), respectively [Ref. l:pp. 150-151). Being a 
recursive procedure, Levinson’s algorithm permits the 
calculation of the (P+l)-st order model parameters by using 
the previously determined P-th order model parameters. The 
matrix form of the algorithm is given by 
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1 




1 




0 


a ( 1 ) 




a (1) 




a (P) 


P+1 




P 




P 


a (2) 




a (2) 




a (P-1) 


P+1 




P 




P 


• 


Z 


• 


- K 


• 


• 




# 

• 


P+1 


# 

• 


a (P) 




a (P) 




a (1) 


P+1 




P 




P 


a (P+1) 




0 




1 


P+1 











(2.64) 



where 



P 

Y a (i)R(P+l-i) 
i = 0 P 

K = . (2.65) 

P+1 P 

y a ( i)R(i) 
i = 0 P 



The P subscript on the "a" parameters specifies the P-th 

order prediction error filter A (z) whereas the index 

P 

identifies the appropriate term in the A (z) polynomial. 

P 

K is called the (P+l)-st order reflection or PARCOR 

P+1 

(partial correlation) coefficient. The PARCOR coefficient 

K represents the true, or direct, correlation between 

P+1 

y(n-P-l) and y(n) with the effects of the intermediate 
variables (i.e., y ( n-P) , y( n-P+ 1 ),..., y( n- 1 ) ) removed. The 
recursive form of Levinson’s algorithm is written as 



a (m)=a (m) -K a (P+l-m) for l<.m_<P, (2.66a) 

P+1 P P+1 P 



and 
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a (P+1) = -K for ra=P+l. (2.66b) 

P+1 P+1 

This shows that the reflection coefficient for each stage P 

is equal to the highest coefficient of A (z). There is a 

P 

one-to-one correspondence between the PARCOR coefficients 

and the coefficients of the transfer function A (z). The 

P 

transfer function or, equivalently, the autocorrelation 

T 

matrix R = E[y(n)y (n)] uniquely determine the reflection 
coefficients [Ref. 7:p. 829]. Taking the z-transform of 
equation (2.66) results in 



A (z) = A (z) - K B (z) (2.67a) 

P+1 P P+1 P 

where 

-1 -2 -P 

A (z) =l+a (l)z +a (2)z + . . . +a (P)z , (2.67b) 

P P P P 

and 

-P -1 

B (z) = z A (z ) , (2.67c) 

P P 

-1 -2 -(P-1) -P 

= a (P)+a (P-l)z +a (P-2)z + . . . +a (l)z +z 

P P P P 

Due to the effective folding about the axis and the shifting 

to the right of the sequence A (z), the polynomial B (z) is 

P P 

called the reverse of A (z). Taking the reverse of both 

P 

sides of equation (2.67) yields 

-1 

B (z) = z B (z) - K A (z) (2.68) 

P+1 P P+1 P 
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Combining equations (2.67) and (2.68) into matrix form gives 







-1 






A ( z ) 




1 -K z 




A ( z ) 


P+1 




P+1 




P 




z 


-1 






B (z) 




-K z 




B (z) 


P+1 




P+1 




P 



Tlie forward prediction error associated with 
predicting y(n) from the previous P samples {y ( n-P) , . . . y (n- 

A 

1)} is written as e (n) = y(n) - y (n) where the subscript P 

P P 

represents the filter order. In z-domain notation, this 
becomes 



E (z) = A (z)Y(z) (2.70) 

P P 

Now, the backward prediction error r (n) is introduced. It 

P 

is defined as the error in predicting ( or actually 
smoothing) y(n-P) from the future P samples {y ( n-P+ 1 , . . . , 
y(n)}. It is written as 



r (n) = y(n-P)+a ( 1 ) y ( n-P+ 1 ) + . . . +a (P)y(n) (2.71a) 

P P P 

or in z-domain notation as 



E (z) = B (z)Y(z) . (2.71b) 

P P 

The optimal backward predictor coefficients minimize the 

2 

mean-squared smoothing error E[r (n)]. Since the optimal 

P 

backward and forward predictor coefficients are mirror 
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2 2 

images of each other, then E[r (n)] = E[e (n)]. That is, 
the backward and forward prediction error vectors have the 
same norms [Ref. 8:p. 101]. 

As will be shown, at each instant n, the backward 
prediction errors are mutually orthogonal (i.e., uncorre- 
lated if assumed to be zero mean) [Ref. l:p. 170]. 
Therefore, it is the backward prediction errors, 

r ( n ) , p=0 , 1 , . . . , M , where M is the filter order, which form 
P 

the new set of basis vectors. To demonstrate that the 
backward prediction errors are mutually orthogonal, let M=3 
and write the corresponding equations for P = 0, 1,2,3 in 

matrix form: 

r (n ) 

0 

r ( n ) 

1 

r (n) 

2 

r (n) 

3 

or 

r^(n) = Ly(n) . (2.72b) 

Note that the first column of L contains the negatives of 
all the reflection coefficients. Now examine the covariance 
matrix of r ( n ) : 



10 0 0 




y (n ) 


a ( 1) 1 0 0 

1 




y(n-i) 


a ( 2 ) a ( 1 ) 1 0 

2 2 




y(n-2) 


a (3) a (2) a (1) 0 
3 3 3 




y(n-3 ) 
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R 



T T 

= E[r(n)r (n)] = E [ Lz( n ) ( L^K n ) ) ] 



(2.73) 



rr 



T T 

= LE[i;(n)x <n)]L 



= LR L 

yy 



Since the normal equations are satisfied within the matrix 

products above, R reduces to a diagonal matrix which 

rr 

verifies that the components of r(n) are uncorrelated. That 
is 



T 

R = LR L = D = diag{E ,E ,E ,E } (2.74) 

rr yy 0 12 3 

where E is the minimum value of the mean-squared prediction 
P 2 

error which is given by E[e (n) •] . It can be shown that 
2 P+1 

E = ( 1-K )E where the recursion is initialized with 
P+1 P+1 P 2 

the value given by E = R (0) = E[y (n)] [Ref. 8:p. 105]. 

0 yy 

Therefore, the prediction error decreases by a factor of 
2 

(1-K ) from one lattice stage to the next. Now, since the 

P+1 

elements of r(n) are uncorrelated, equation (2.72) is 

equivalent to the Gram-Schmidt orthogonalization of the data 

vector y(n). Also note that rewriting equation (2.74) as 
-1 -T 

R = L DL corresponds to a LU-Cholesky factorization of 

yy 

the covariance matrix of y(n). 

To complete the derivation of the lattice recursion 
equations, multiply both sides of equation (2.69) by Y(z) to 
get 
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(2.75) 







-1 






E (z) 




1 -K z 




E (z) 


P+1 




P+1 




P 






-1 






E (z) 




-K z 




E (z) 


P+1 




P+1 




P 



These equations are transformed into -the time domain to 
obtain the final result: 





(n) = e (n) - 


K r .( n- 1 ) 


( 2 . 76a) 


P+1 


P 


P+1 P 






(n) = r (n-1) 


- K e (n) . 


( 2 . 76b) 


P+1 


P 


P+1 P 





These equations, which are recursive in both time and order, 
define the signal flow graph of the analysis lattice filter, 
shown in Figure 2.6. For a given time instant n, the equa- 
tions are evaluated recursively in order. The inputs into 

the first stage of the lattice filter are e (n) = r (n) = 

0 0 

y(n). The successive stages of the filter develop the 

successively higher order forward and backward prediction 

errors. The output from the final stage yields the desired 

M-th order forward prediction error, while all lower order 

prediction errors are available at intermediate stages. 

Since the backward errors are generated from the y data 

vector, the analysis lattice filter actually implements 

the Gram-Schmidt orthogonalization ; all that is required to 

implement the lattice are the reflection factors 

{K ,K , . . . ,K } . 

12 M 
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Figure 2.6 Analysis Lattice Filter 
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Equation (2.76) can be manipulated to obtain the 

lattice form of the synthesis filter H(z) = 1/A (z): 

M 

e (n) = e (n) + K r (n-1) (2.77a) 

P P+1 P+1 P 

r (n) = r (n-1) - K e (n) . (2.77b) 

P+1 P P+1 P 

Figure 2.7 shows the corresponding signal flow graph. When 

the input to the synthesis filter is the forward prediction 

error sequence e (n), the output is the original sequence 

P 

y(n). In order for the synthesis filter to be stable and 

causal, all M zeros of the prediction-error filter A (z) 

M 

must lie within the unit circle in the complex z-plane. A 

necessary and sufficient condition for all the zeros to be 

inside the unit circle is that the magnitude of each of the 

reflection coefficients {K ,K } be less than one. 

12 M 

[Ref. l:pp. 168-169] 

The reflection factors can be evaluated by several 
methods. The various methods arise due to different 
definitions of the optimality criterion. The criterion 
used here of minimizing the mean-squared prediction errors 
leads to what MAKHOUL calls the forward and backward methods 
[Ref. 9]. They are, respectively: 

2 

K = E[e (n)r (n-1)] / E[r (n-1)] (2.78a) 

P+l,e P P P 

2 

K = E[e (n)r (n-1)] / E[e (n)] . (2.78b) 

P+l,r P P P 
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Figure 2.7 Synthesis Lattice Filter 



54 





2 2 

Assuming stationarity , and since E[r (n) ] = E[e (n) ] as 

P P 

previously discussed, then the forward and backward 

reflection coefficients are equal. That is K = K = 

P+1 P+l,e 

K . The Schwarz inequality implies that [K ! ^ 1 for 

P+l,r P 

each P=1,2,...,M [Ref. 8:p. 104]. An alternate technique for 

calculating the reflection coefficients is the geometric 

mean method. It was introduced by ITAKURA and SAITO in 

their work of developing a digital filter structure for 

time-domain speech analysis [Ref. 10]. The corresponding 

PARCOR coefficient is given by 



E[e (n)r (n-1)] 

P P 

K = (2.79) 

P+1 2 2 1/2 

[E[e (n)]E[r (n-1)]} 

P P 

These PARCOR coefficients are guarenteed to have magnitudes 

less than one [Ref. 7:p. 840]. 

A third method was used by BURG in the maximum- 

entropy method of spectral estimation [Ref. 11]. This is 

the harmonic-mean method. The harmonic-mean method seeks to 

minimize the sum of the forward and backward prediction 

2 2 

error variances, E[e (n)] + E[r (n-1)]. Minimizing this 

P P 

sum with respect to the reflection coefficients results in 



55 



W t HKO Ci (J fc. U M f O vj V t:. f< l>l e. f>( I t. r c. im oc. 



(2.80) 



2E[e (n)r (n-1)] 

P P 

K = 

P+1 2 2 

E[e (n) + E[r (n-1 ) ] 

P p 

Again, the Schwarz inequality verifies that K has magnitude 

P 

less than one, guaranteeing that the synthesis filter is 
causal and stable [Ref. l:p. 189]. 

4 . The Generalized (Analysis) Lattice Filter 

The preceding development of the analysis lattice 
filter equations assumed that the data sequence {y(n)) 
represented a time sequence with stationary statistics. In 
this section, a more general linear prediction problem 
will be considered. No special assumptions are made 
concerning the data. The data values need not be delayed 
versions of each other; they need not even represent a time 
sequence. This is the approach taken by LENK in developing 
the generalized order update equations [Ref. 12:p. 85]. The 
resulting generalized form of the lattice filter makes it 
suitable for multidimensional and nonlinear signal 
processing applications. 

Definitions associated with the normalized form of 
the generalized lattice filter will now be introduced. 
First, the components of the length (M+1) column vector 
[y'^ ], representing a single realization of the random 
process Y, are designated by y^ , A= 0,1,..., M. The forward 
prediction error in estimating the (n+l)-st element from 
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the previous N elements of the data vector is written as 
N M N X 

e = Z ^ (2.81) 

n+1 X-0 ^ 

where the length (M+1) row vector of coefficients is 
given by 



N N N N 

fh ( n+ 1 ) ] = [ 0 , . . . , 0 , -h ,-h , . . . , -h , 1 , 0 , . . . , 0 ] . ( 2 . 82 ) 

^ n-N+1 n-N+2 n 

The backwards prediction error associated with predicting 
n-N 

y from the next N elements of the data vector is given by 



N 



M . N 



= V h ( n-N ) y 

X 



(2.83) 



n-N A=0 

where the associated coefficients are given by the vector 



N 



N 



N 



N 



[h ( n-N ) ] = [ 0 , . . . 0 , 1 , -h ,-h , . . . , -h , 0 , . . . , 0 ] . ( 2 . 84 ) 

^ n-N+1 n-N+2 n 



The norm of the forward prediction error is defined as 
N N 2 1/2 

Me M = [E{(e ) } ] . (2.85) 

n+1 n+1 

The norm of the backward prediction error is defined in a 

similar manner. Now, the normalized forward and backward 

N-th order prediction errors are defined by 

N N N M N 

e = e / Me M = ^ a (n+l)y , (2.86a) 

n+1 n+1 n+1 ^ 

and 

_N N N M N 

r = r / Mr “ Z t>s(n-N)y (2.86b) 

n-N n-N n-N \=0 ^ 
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where 



the normalized prediction error coefficients are 



defined as 

N N N 

a (n+1) = h (n+1) /lie | | (2.87a) 

n+1 



and 

N N 

b (n-N) = h (n-N) / ||r || . (2.87b) 

A A 

[Ref. 12:pp. 85-87] 

Using the normalized form of the generalized 
Levinson algorithm, LENK demonstrated that equations (2.87) 
could be updated recursively through the relation 



N+1 




N 


a . ( n+ 1 ) 


n 


a ( n+ 1 ) 


X 


= 0(K > 


X 


N + 1 


N+1 


N 


b ^ (n-N) 




b- (n-N) 



( 2 . 88 ) 



where the partial correlation ( PARCOR ) coefficient is 
n N N 

K = E{e r } , (2.89) 

N+1 n+1 n-N 



and where 



n 



0(K ) 

N+1 



n 2 

1 - (K ) 
N+1 



- K 



n 



- K 



n 

N+1 

1 



N+1 



(2.90) 



The recursion is started for order N=0 with the initial pre- 
diction error values for A=0,1,...,M given by the vectors 
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0 n+1 

[a^(n+D] = [0,...,0, l/||y || ,0,...,0] (2.91a) 

0 n 

[b^(n)] = [0,...,0, l/||y n, 0,...,0] . (2.91b) 



Multiplying both sides of equation (2.90) by the data vector 
[y'^ ] yields the desired error order update equations: 



_N+1 

e 

n+1 

_N+1 

r 

n-N 



e<K " ) 

N+1 



_N 

e 

n+ 1 

_N 

r 

n-N 



(2.92) 



The corresponding signal flow graph for a single lattice 
filter section is shown in Figure 2.8. Figure 2.9 depicts a 
third order generalized lattice filter. [Ref. 12:pp. 94-99] 
LENK then proved that the reflection coefficients 
(i.e., PARCOR coefficients) could be calculated directly 
from the data vector’s correlation matrix by utilizing the 
generalized Schur algorithm. In LENK’s notation, the 
components of the correlation matrix are written as R = 
E(y^ y^ }. Then the reflection coefficients are given by 



n 

K = 
N+1 



n-N 

a (n+1) 
N 



P 



n-N 

(n-N) 

N 



(2.93 ) 



where 
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N+1 
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a. Single Stage of Generalized Lattice Filter 



N+1 




N+1 

e 

n+1 



b. Alternate Representation 



Figure 2.8 Generalized Lattice Filter 
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Figure 2.9 3— rd Order Generalized Lattice Filter 
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rttrrti^uuK-iLU 



5 




< 








D 

X, 



and 



X N X M N 

ly (n+l ) = E{e y } = y a (n+1 )R 
N+ 1 n+l |4 =0 




(n-N) 



NX M N uA 

E{” y } = y b (n-N)R 
n-N =0 



( 2.94a) 



(2.94b) 



Equations (2.94) can be updated through the recursion 
relation 



A 




A 


y (n+l) 




y (n+l) 


N+l 


n 


N 




= (K ) 




^ A 


N+l 




n (n-N) 

^N+1 




n (n-N) 
N 



(2.95) 



To start the recursion at order N=0, the values of the 

parameters for ^=0,1,...,M are given by the vectors 

\ n+l -1 (n+l)0 (n+l)l (n+l)M 

tcy (n+l) ] = | ly II [R ,R ,...,R ] (2.96a) 

0 

X n-1 nO nl nM 

[/? <n)] = lly M [R ,R ,...,R ] . (2.96b) 

•^0 



[Ref. 12:pp. 100-101] 

5 . Lattice Filter Advantages 

The lattice filter has several advantages compared 

to the direct, or transversal, form of the prediction error 

filter A (z). First of all, due to the built-in 
M 

orthogonali zation incorporated into the lattice, successive 
lattice stages are decoupled. Therefore, the reflection 
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coefficients K , P=1,2,...,M are independent of the filter’s 
P 

final order. The M-th order least squares prediction filter 

contains all the prediction error filters of lower order. 

Restated, the first P sections of the M-th order filter form 

the P-th order prediction filter. Therefore, lattice stages 

may be added or subtracted from the existing lattice filter 

without having to recalculate the already determined 

reflection coefficients. In contrast, when the order of the 

transversal filter is changed, all the A (z) prediction 

M 

error filter coefficients must be recalculated. As a 

consequence, to specify all filter orders up to M requires 
storage of M(M+l)/2 predictor coefficients, while only M 
reflection coefficients need to be stored. [Ref. 7:p. 830] 

Another advantage of the lattice over the 

transversal filter is that the output of each lattice stage 
is a least-squares prediction error. Therefore, if the 
desired order of the predictor is not known in advance, the 
output of the various stages can be monitored to determine 
what filter order is adequate. [Ref. 8:p. 106] 

Due to the quantization inherent in implementing 
digital filters, round off noise is introduced.- Studies 
have shown that the performance of the lattice filter is 
far superior to that of transversl filters for finite word 
length computations. Furthermore, the filter zeros are less 
sensitive to quantization errors in the reflection 
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coefficients than for quantization errors in the transversal 
filter coefficients. Finally, since the magnitude of the 
reflection coefficients is less than one, it simplifies the 
task of establishing the quantizer overload point. [Kef. 

8:p. 106] 

Another advantage the lattice filter holds over the 
transveral filter is that the lattice filter is minimum 
phase (i.e., all its zeros are inside the unit circle so 
that the synthesis model is stable) if and only if the 
magnitude of the reflection coefficients is less than one. 
There is no comparable test for the transversal filter 
coefficients to determine whether the synthesis model is 
stable. [Ref. 8:p. 106] 

6 . Linear Lattice Filter Applied to Deconvolution 

The transfer function of the lattice filter is 

determined by the reflection coefficients. In turn, the 

values of these reflection coefficients are uniquely 

determined by the prediction error filter transfer function 

A (z), or equivalently by the autocorrelation sequence R(k), 
M 

k = 0,1, ....M [Ref. 7:p. 829]. Similarly, equations (2.70) 

and (2.71b) indicate that the transfer -functions from the 

input y(n) to the outputs e (n) and r (n) are A (z) and 

M M M 

B (z), respectively. Therefore, the analysis lattice filter 
M 

is equivalent to the whitening filter A (z); that is, it is 

M 

the inverse of the system model H(z) = 1/A (z). This is the 

M 
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basis for using the lattice filter in a deconvolution 
application . 

In order to recover the input signal x(n) from the 

system’s output signal y(n), it is necessary to determine 

the inverse of the system transfer function H(z). The 

initial step is to conduct an "identification experiment" to 

estimate the system’s parameters (either the autoregressive 

filter coefficients or the lattice filter reflection 

coefficients). To accurately estimate the parameters, the 

chosen experimental input signal must be sufficiently rich 

in frequency content, or more formally, persistently 

exciting so that it excites all the modes of the system. A 

sequence is said to be persistently exciting of order n if 

its (n X n) autocorrelation matrix is nonsingular [Ref. 

13:pp. 70-71]. By using the techniques presented in section 

D.3 of this chapter, the resulting output sequence y(n) can 

be processed to yield the desired reflection coefficients. 

When the analysis lattice filter is implemented with these 

coefficients embedded in its structure, the lattice becomes 

the inverse of A (z) = 1/H(z). This procedure is depicted 

M 

in Figure 2.10. 

The application of linear lattice filters to decon- 
volution was simulated using the FORTRAN programs LININV, 
ATOCOR , LEVIN, AND LATICE. These programs are provided in 
Appendix A. The simulation was conducted for a second 
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noise 




noise 



Figure 2.10 Deconvolution Using Lattice Filter 
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order, autoregressive, LTI system defined by the equation 
y(n) = x(n) - (0.6y(n-l) + 0.08y(n-2)) . As previously 
mentioned, modeling the system was the first step in the 
deconvolution process. In order to identify the "unknown" 
parameters, the system was excited by a zero mean, unity 
variance, Gaussian white noise sequence. The output 
sequence y(n) was processed by subroutine ATOCOR to 
calculate the components of the autocorrelation matrix R 

yy 

Then subroutine LEVIN implemented Levinson’s algorithm to 
evaluate both the prediction error filter coefficients 
and the associated reflection coefficients from the autocor- 
relation matrix. The actual output of subroutine LEVIN is 
the lower triangular L matrix described in equation (2.72); 
the i-th row of L contains the i-th order prediction error 
filter coefficients listed in reverse order, and the first 
column contains the negative of the reflecti'on coefficients. 
Once the reflection coefficients corresponding to 1/H(z) 
were calculated, they were embedded in subroutine LATICE 
which implemented the analysis lattice filter equations 
(2.76). With the inverse filter of H(z) now available 
(i.e., the analysis lattice filter), H(z) was driven by an 
"unknown" signal x(n). The output of H(z) was then fed into 
the lattice filter. An approximation to x(n) was recovered 
at the forward error output signal from the last stage of 
the lattice. 
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Simulations were run both with and without the 



presence of measurement noise. Table 2.1 displays the 
resulting L matrices for one simulation. The case of no 
measurement noise is shown in Figures 2.11 through 2.14. 
Figure 2.11 is a plot of the input signal x(n), Figure 2.12 
depicts the system output y(n), and Figure 2.13 shows the 
lattice filter output x(n). As can be seen, the results of 
the inverse filtering were excel lent--the x and x curves are 
identical. This is verified by Figure 2.14 which shows that 
the mean-square error between x(n) and x(n) is nearly zero. 
The running average mean-square error was calculated using 
the equation 



MSE(n) 



n 

(1/n) y 

i = l 



(x(i) - 



x(i 



)) 



(2.97) 



The simulation was then repeated for a nonzero 
measurement noise. Here, the added measurement noise was a 
zero mean, 0.0025 variance, Gaussian white noise sequence. 
Using the same input as for the previous case, the outputs 
of the system and lattice filter are shown in Figures 2.15 
and 2.16, respectively. Figure 2.17 is a plot of the mean- 
square error between x(n) and ^(n). The results in Table 2.1 
show that even with the presence of measurement noise, the 
system parameters were still accurately identified. The 
lattice filter recovered the basic shape of x(n), however, 
it could not remove the additive measurement noise. 
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the output of the lattice filter is a noisy 



Therefore , 

or "fuzzy" version of x(n). Additional 
to remove the noise component of 
presented here are representative of 
simulations involving other stable, 
systems . 



filtering is required 

A 

x(n). The results 
those obtained for 
autoregressive, LTI 
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TABLE 2 . 1 



RESULTS OF MODELING A LINEAR SYSTEM 



Modeling Problem : 

-1 -2 

System: H{z) = 1/[1 + 0.6z + .08z ] 

System input: Gaussian white noise, N(0,1) 
Number of white noise realizations used: 25 

Number of points per realization: 5,000 



Modeling Results : 



a. No Measurement Noise 



L = 



1 0 
0.555 1.0 

0.077 0.598 



0 

0 

1 



Order , P Ep 



Er 



0 1.445 

1 0.999 -0.555 

2 0.993 -0.077 



b. Measurement Noise is N(0, 0.0025) 



L 



1 0 

0.555 1.0 

0.077 0.598 



0 

0 

1 



Order , P Ep Kp 

0 1.449 

1 1.002 -0.555 

2 0.996 -0.077 
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Figure 2.11 System Input, x(n) 
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Figure 2.12 System Output, y(n) 
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Figure 2.13 Lattice Filter Output, x(n) 
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Mean-Square Error Between x(n) and x(n) 
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Figure 2.14 
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Figure 2.15 System Output Plus Measurement Noise 
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Figure 2.16 



Lattice Filter Output, x(n) 
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Figure 2.17 Mean-Square Error Between x(n> and x(n) 
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A. INTRODUCTION TO NONLINEAR SYSTEMS 

While the previous chapter dealt solely with linear 
systems or plants, this chapter will address the inverse 
filtering problem involving nonlinear systems. The chapter 
starts with a brief introduction to modeling nonlinear 
systems. Then it quickly proceeds to extend the generalized 
linear lattice filter results of section II. D. 4 to a 
nonlinear analysis lattice filter. The nonlinear lattice 
filter is discussed in detail. To conclude, results from 
numerous simulations involving the lattice in deconvolution 
applications are presented. 

System linearity is defined in terms of the principle of 
superposition. If the rule by which the system transforms 
the input x(k) into the output y(k) is represented by the 
operator T, then the system is said to be linear if and only 
if 

T[ax (k) + bx (k)] = aT[x (k)] + bT[x (k)] (3.1) 

12 1 2 

= ay (k) + by (k) 

1 2 

for arbitrary constants a and b. The system is nonlinear 
if equation (3.1) is not satisfied. Estimation of the 
parameters of a nonlinear system is a complex problem. 
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One approach to the parameter estimation problem for 

nonlinear filters is the Bayesian approach. This approach 
leads to various approximate solutions for the parameter 

estimates. In the Bayesian approach, the parameters are 
considered to be random variables. The parameter vector h 
is considered a random vector with a probability density 
function p(h). Introducing the observation vector y(t) and 
the input vector x(t), both of which contain data up to 
time t, the a posteriori probability density function for h 
is p ( h I t ) , X ( t ) ) . One possible choice for the estimate of 

A 

the h vector is to select the conditional mean h ( t ) = 

E [h I y,( t ) ^( t ) ] . Selecting this as the estimate minimizes the 

variance of the parameter estimation error. Another 

A 

possible choice is to select the h(t) which maximizes 
p( h I y;( t ) , X ( t ) ) . This most likely value is known as the 
maximum a posteriori (MAP) estimate. Finally, the Bayesian 
approach also leads to the extended Kalman filter for 
nonlinear state estimation problems. [Ref. 13:pp. 32-41] 

Another popular method for modeling nonlinear systems is 
based on the Volterra series. The series is named after the 
mathematician Vito VOLTERRA. The first person to apply the 
series to nonlinear systems was Norbert WIENER. If the 
"black box" approach is taken towards a nonlinear, time- 
invariant system, the relationship between the input x(t) 
and the output y(t) can be represented by the Volterra 
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series : 



y(t) 



-oo 

'oo 



= f h (T 

■i 

iJL 



)x(t-T )dT 
1 1 



+11 h (T ,T )x(t-T )x(t-T )dT dT 

.CO 2 1 2 1 2 12 

'00 

+111 h (T ,T ,T )x(t-T )x(t-T )x(t-T )dT dT dT 
.oo 3123 1 2 3123 



^oo roo 
f / . . ./ h (T 

Ico 1 



, ...,T )x(t-T )...x(t,T )dT . . .dT 
n 1 n 1 n 



+ . . . (3.2) 

where n = 1,2,... and h (T ,...,T ) =0 for any T < 0, 

n 1 n j 

j = l,2,...,n . The functions h (T , . . . ,T ) are called 

n 1 n 

Volterra kernels of the system. This equation is a 

functional series. That is, it performs an operation on the 

function x(t) which results in a number for y(t). If the 

n-th order Volterra operator, H , is introduced where 

n 



y (t) = H [x(t)] (3.3) 

n n 

(T ,...,T )x(t-T )...x(t-T ) dT ...dT , 
n 1 n 1 n 1 n 

then the series can be expressed in operator notation as 




y(t) = H [x(t)] + H [x(t)J + ... + H [x(t)l +... (3.4) 

1 2 n 

00 oo 

= y H[x(t)]= y y(t) . 

n=l n n=l n 
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Figure 3.1 is a graphic representaion of this equation. 

[Ref. 14:pp. 7-9] The H operator is a linear operator, 

1 

H is a quadratic operator; and, in general, the H operator 
2 n 

involves the terra x(t) to the n-th power. 

B. GENERALIZED NONLINEAR LATTICE FILTER 

1 . Introduction 

In this section, it will be demonstrated how the 
generalized lattice filter introduced in section II. D. 4 can 
be applied to modeling nonlinear systems. The development 
will start by looking at the discrete form of the Volterra 
series. Based upon this series an alternate tensor notation 
representation will be introduced. The results of section 
II. D. 4 will then be extended to handle a two-dimensional 
field of data which will result in the generalized nonlinear 
lattice filter. 

2 . Nonlinear Lattice Filter Development 

The discrete form of the Volterra series of equation 
(3.2) is given by 

y(k) = h + h (n )x(k-n ) (3.5) 

0 11 1 

+ h (n ,n )x(k-n )x(k-n ) + ... 

2 12 1 2 

Using LENK’s tensor notation, an equivalent form of equation 
(3.5) is given by 
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y (t) 










Figure 3.1 Volterra Series Representation 
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y(k) 



• • ♦ 



, where 



(3.6a) 



y(k) = 
x^\ k ) 


A 1 

H+Hx +H XX + ... , where (3.6a) 

\ 

T 

= [ X ( k ), X ( k- 1 ),..., X ( k- Xj ),..., X ( k-N ) ] (3.6b) 



for A; = 0,1,..., N. As an example, if N=1 and equation 
(3.6a) is truncated to the three terms shown above, then 
this equation can be rewritten as 



y( k ) = 


H + H x(k) + H x(k-l) + H x(k)x(k) 
0 1 00 

+H x(k)x(k-l) + H x(k-l)x(k) 

01 10 




+ H x(k-l)x(k-l) . (3.7) 

11 



Based on this tensor notation, LENK introduced an alternate 
representation for the nonlinear system. Instead of 

defining the components of the vector x (k) as the present 
and past values of the signal x(k) as in equation (3.6b), 
the components are redefined in terras of x(k) raised to the 
power for A*, = 0,1,2,.. .,N. That is 



Ai , , V 

X (k) 


0 12 NT 

= [x (k),x (k),x (k),...,x (k)] (3.8) 

2 NT 

= [ 1 ,x(k) ,x (k) , . . . ,x (k) ] 



for A'l = 0,1,2,...,N. Now the output of the nonlinear system 
is defined by 

y(k) = x"^(k) . . .x^'(k-N)H . . . .^ (3.9) 

Ao Ari 
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I 



!* 

1^') 







':/C 



o 

-'i j 







'■i 



for^Q,A|)... = 0,1,2, ...,P. The variable P gives the order 

of the filter and N defines the filter’s finite memory. The 

term H ^ • • --s plays a role similar to that of the Volterra 
Ao Ah 

kernel; it can be considered a (P+1) -order tensor. To 
clarifly the meaning of this notation, the following example 
with N=1 and P=2 is provided: 



y(k) = lk)x '(k-1) . 

2 

= H +H x(k)+H x(k-l)+H x(k)x(k-l)+H x (k) 
00 10 01 11 20 
2 2 2 
+ H X (k-l)+H X (k)x(k-l)+H x(k)x (k-1) 

02 21 12 



2 2 

+ H X (k)x (k-1 ) . (3.10) 

22 

[Ref. 12:pp. 39-48] 

The above nonlinear model for y(k) is a moving 
average (MA) model: the output is defined in terms of the 

input signal. The next step is to establish an 
autoregressive (AR) model where y(k) is estimated in terms 
of its past values. This type of model is useful when the 
input signal is not readily available. An autoregressive 
model is obtained by redefining the observation vector in 
terms of y(k) vice x(k). Then the AR model is given by 

y(k) = y (k-l)...y (k-N)H ..., (3.11) 

A, An 

for = 1,2,...,P. If the system is driven by a white 

noise signal u(k), and if the system’s parameters are 
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A 

approximated by H . .*. , then the estimation error is 

Ac 

e(k)=y(k)-y(k) = [y^\k-l) . . .y‘^1k-N)H ... + u(k)] 

A\ An 

- [ y'^\k-l) . . .y^'(k-N)H, . . .;l ] (3.12) 

Aj 



If the system parameters are known exactly, then the estima- 
tion error and the white noise sequences are equal — that is 
e(k)=u(k). One note concerning the nonlinear AR model: 



unlike the linear AR model whose stability is easily 



determined 


( i .e. , 


the system is 


stable 


if 


all its 


poles 


lie 


within the 


unit 


circle in the 


complex 


z-plane ) , it 


is 


difficult 


to judge for which 


class 


of 


inputs 


the 


AR 


nonlinear 


system ’ 


s output will 


remain 


bounded . 


This 


is 



because the order of the nonlinearity increases with time. 
[Ref. 12;pp. 68-69] 

Using this alternate tensor form of the AR nonlinear 
system, along with the generalized lattice of section 
II. D. 4, the nonlinear lattice filter will be developed. To 
simplify the discussion, the filter’s finite memory will be 
restricted to N = 2. Then the product Y(k) = [ y^'( k- 1 ) y^^ k-2 ) ] 
forA\,Aj_ =0,1,..., P forms a second order tensor, or a two- 
dimensional data field given by: 
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Y(k) 



P 

1 y(k-2) ... y (k-2) 

P 

y(k-l) y(k-l)y(k-2) ... y(k-l)y (k-2) 

2 2 2 P 

y (k-1) y (k-1 )y(k-2) . . . y (k-l)y (k-2) 



(3.13) 



P P P P 
y (k-1) y ( k-1 ) y ( k-2 ) . . . y (k-l)y (k-2) 

The types of systems that this nonlinear lattice 
structure can model exactly are of the form shown in Figure 
3.2 . The nonlinear combinations block forms a weighted sum 
of the cross-products and powers of the input variables 
y(k-i), for i=l,2,...,N. As an example, with N=2 and P=2 the 
general equation for the system’s output when excited by the 
input x(k) is given by: 



y(k) = x(k) - { H + H y(k-l) + H y(k-2) 

11 21 12 
2 2 

+ H y(k-l)y(k-2) + H y (k-1) + H y (k-l)y(k-2) 

22 31 32 

2 2 
+ H y (k-2) + H y(k-l)y (k-2) 

13 23 

2 2 

+ H y (k-l)y (k-2) } (3.14) 

33 



It is also assumed that the system is time-invariant; 

otherwise the model changes with time k. In order to define 

2 

the forward and backward prediction errors, the (P+1) 
elements of the two-dimensional data field must by converted 
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Figure 3.2 Nonlinear System Model 
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into a one-dimensional sequence. The ordering chosen by 
LENK is listed below: 

order = {( P , P ),( P-1 , P) ,( P, P-1 0 , P) ,< P , 0 ) , 

(P-1,P-1) , (P-2,P-1) , (P-l,P-2) , . . . , (0,P-1) , 
(P-1,0),. ..,a,l),(0,l),(l,0), (0,0)} , (3.15) 



where the components of the data field are identified by the 

indices (m,n) for m,n = 0,1,... ,P. The elements of the 

2 

ordered set are numbered consecutively from 0 to (P+1) -1. 

The notation (m,n)-q is used to identify the q-th element 

prior to the element (m,n) as referenced to the ordered 

sequence. This notation will be further abbreviated to 

mn-q. Now the (q-l)-order, normalized, forward error 

associated with predicting the value of the element 
m n 

y (k-l)y (k-2) from the preceding (q-1) elements is given by 



A, Ar 

e = a (m,n)y (k-l)y (k-2) 
mn A( Aj. 



(3.16) 
q-1 



forA»,Aj =0,1,2, ...,P. Note that the coefficients a^ -v can 

Aa, 

be thought of as the components of a second order tensor. 

q-1 

The coefficient a. ^ is equal to zero when the indices 
(A| »^i) do not correspond to the (q-1) elements preceding 
(m,n) in the ordered sequence (i.e.f when (Aj,Aj,) ^ (ni,n) 
or when (AviAj^) <. mn-q ). Also, when - (m,n), then 



q-1 q-1 

a (m,n) = 1 / I I e | | 

mn mn 



(3.17) 
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Similarly, the normalized backward prediction error in 
estimating y(mn-q) from the next (q-1) points is given by 



_q-i q-1 X, Az 

r ~ (®n-q)y (k-l)y (k-2) 

mn-q 



(3.18) 



q-1 



for A| , Xx =0 » 1 » 2 , . . . ,P . The coefficients b. ,(mn-q) equal 

X\Ax 

zero when the indices do not correspond to the (q-1) 

elements following (m,n) (that is, when (Ai»A 2 _) < (mn-q) or 
(A|,A 2 .) > (m>n) ). When (Ai»A 2 ,) = (m,n), then 



q-1 

b (mn-q) = 
mn-q 



q-1 

1/ I Ir II. 
mn-q 



(3.19) 



[Ref. 12:pp. 145-146] 

Using the normalized nonlinear Levinson algorithm, 
LENK shows that the nonlinear prediction errors can be 
updated in order through the recursion relation 







^q-1 


e 




e 


mn 


mn 


mn 




II 






q 


_q-l 


r 




r 


mn-q 




mn-q 



(3.20) 



where 



mn 



0(K ) = — 

* ^ 



mn 2 
(K ) 

q 



mn 



-K 



- K 



mn 

q 

1 



(3.21) 
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and the unique reflection coefficients are given by 



mn 



q-l q-1 



K 



E{e 



r } 
mn mn-q 



(3.22) 



q 



[Ref. 12:pp. 146-147] 



A deficiency remains to be corrected: the goal is 



to estimate y(k), but there is no y(k) term in the two- 
dimensional data matrix of equation (3.13). This problem is 
solved by adding another channel to the lattice structure 
and by exploiting the orthogonality of the backward pre- 
diction errors. Since the backward prediction errors leaving 
the top row of the lattice filter (Figure 2.9) are uncor- 
related, they can be used in a Fourier series to estimate 

y(k). These backward prediction errors can be formed into a 

2 

length L = (P+1) vector defined as 



where the subscript is in the form mn-q. Now, the error in 
estimating y(k) using the data in the Y(k) tensor is 
calculated from 



A" 

( m , n ) - x' 



0 12 
] = [r ,r ,r 

00 00-1 00-2 



, . . . , r 



L-1 T 

] (3.23) 

00-L+l 



e 



k 



L 



y(k) 




(3.24) 



where the Fourier coefficients, K / , are given by 
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} ] . (3.25) 



0 1 L-1 

[K,]=[E{y(k)r },E{y(k)7- E { y( k )7 

A 00 00-1 00-L+l 

The resulting generalized nonlinear analysis lattice 

structure is depicted in Figure 3.3. [Ref. 12:pp. 147-150] 

The nonlinear lattice structure will now be examined 

more closely. As can be seen in Figure 3.3 , the inputs to 

the analysis filter are the normalized values of y(k) and 
2 

the (P+1) +1 ordered components of the two-dimensional data 

field (equation (3.13)). The ordering of these inputs is in 

accordance with equation (3.15). At a given time k, these 
2 

(P+1) +2 inputs are evaluated and inserted into the left 
side of the lattice structure. The lattice calculations are 
conducted by starting at the top row and moving downward 
along the northwest-southeast diagonal, and then advancing 
to the top of the next diagonal and so on. When all the 
lattice calculations for time k have been completed, the 
process is repeated for time k+1. Just as in the case of 
the linear lattice filter, the PARCOR coefficients at each 
lattice section act to decorrelate the two input signals. 
The backward error signals exit at the top of the lattice 
structure, and the forward error signals exit at the right 
side. The output of interest for inverse filtering is the 
forward error signal from the top row of the filter (i.e., 
the row corresponding to the y(k) input signal). This is 
the error arising from the estimation of y(k). 
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Figure 3.3 Quadratic Nonlinear Lattice Filter 
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3 . The Nonlinear Lattice Applied to Deconvolution 

The generalized nonlinear lattice filter can be 
applied to the inverse filtering problem just as the linear 
lattice filter was in section II. D. 6. The first problem is 
one of system identification and parameter estimation. Once 
the model parameters have been estimated, they are embedded 
in the nonlinear filter lattice structure. Then, the 
nonlinear lattice represents the inverse of the system’s 
transfer function. As will be shown, the nonlinear lattice 
filter can model both linear and nonlinear systems-- the 
linear system is just a special case of the nonlinear 
system. This inverse filtering algorithm was implemented by 
the FORTRAN programs NLMAIN, NLCLAT, SCHUR, NORMS, NLLAT, 
and URAND. These programs are listed in Appendix B. Now, 
this inverse filtering procedure will be described in more 
detail . 

In order to identify the model parameters (i.e., the 
PARCOR coefficients), the system was excited with zero mean, 
unit variance white noise sequences. Both Guassian and 
uniform noise distributions were used with good results. 
One difficulty encountered in generating the nonlinear data 
was ensuring that the output of the postulated nonlinear 
system remained bounded for the input noise signal. For the 
simulations, the filter and system were constrained to a 
finite memory of at most two delays (N=2), while the 
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nonlinear order was allowed to vary from P=0'to P=4. (For N 

greater than two, the problem of ordering the elements of 

the N-dimensional Y tensor increases in complexity. ) The 

system output sequence y(k) was then processed by 

subroutines NLCLAT and SCHUR, which determined the 

corresponding autocorrelation matrix and the partial 

correlation coefficients, respectively. In an effort to 

improve the accuracy of these calculations, noise sequences 

of up to 5,000 points were used. Additionally, the PARCOR 

coefficients were averaged over as many as 50 realizations 

of the input noise random process. Through trial and error, 

it was found that the best inverse filtering results were 

obtained when the resulting reflection coefficients were 

truncated to two decimal places. The output of subroutine 

2 2 

SCHUR is a ((P+1) +1 x (P+1) +1) upper triangular matrix of 
reflection coefficients. This matrix is simply overlaid 
atop the upper triangular shaped nonlinear lattice structure 
to place the reflection coefficients at the correct filter 
sections . 

With the reflection coefficients calculated and 

embedded in the filter, we were able to use the lattice in 

an inverse filtering application. To recover the input 

signal x(k), the system’s output signal y(k) was processed 

by subroutines NORMS and NLLAT. The function of NORMS was 

2 

to normalize the (P+1) +1 input signals into the lattice. 
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Then subroutine NLLAT, using the previously calculated 
reflection coefficients, implemented the lattice structure 
and carried out the lattice filter calculations. The forv?ard 
error signal out of the top row of the lattice yielded the 
normalized estimate of x(k). Since the inputs to the lat- 
tice filter are all normalized, the outputs must be 
denorraal i zed . Therefore, since the input into the top row of 
the lattice is divided by the norm of y(k), the forward 
error signal from the top row of the lattice is multiplied 
by the norm of y(k). (Note that if y(k) is a zero mean 
sequence, then this norm is equivalent to its standard 
deviation.) Also, it was found that to achieve a good 
estimate, it was necessary to further scale this error 
signal by dividing it by the norm of the noise generated 
sequence y ( k ) . 

In order to verify the accuracy of the reflection 
coefficients, the noise generated output signal, y(k), was 
passed through the inverse lattice filter to see how well 
the filter whitened this sequence. The effectiveness of the 
whitening filter was evaluated by examining the mean-square 
error (MSE) between the white noise input signal, x(k), and 

A 

the lattice filter's output, x(k). The running average 
mean-square error was calculated using the equation 



MSE(k) 





(x(i) 



A 

X(i 



)) 



( 3 . 26 ) 
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J 







I;, I 



!' 1 





Since the noise input, x(k), has unit variance, the MSE 

A 

plot essentially provides a percentage error between x(k) 
and x(k). Plots of the MSE are included in the 
simulation results. Having demonstrated that the lattice 
filter represented a good inverse of the system H(z), the 
system was then driven by a known input signal x(k). To 
recover an approximation to this signal, the corresponding 
system output was passed through the lattice filter, and the 
resulting lattice filter output was denormalized and 
rescaled as previously discussed. Simulation results for 
various systems are presented in the following section. 

4 . Inverse Filtering Simulation Results 

In this section, the previous modeling and inverse 
filtering procedures are implemented and applied to various 
linear and nonlinear systems. Unless otherwise noted, the 
reflection coefficients for each of the systems were 
determined by using twenty-five realizations (5,000 points 
each ) of the zero mean, unit variance white noise random 
process as the input excitation signal. As previously 
mentioned, the mean-square error between this white noise 
input and the lattice filter’s output was plotted to 
evaluate the lattice filter’s inverse filtering performance. 
After the system was modeled, it was driven by the signal 
x(k) which consisted of ramps, pulses, and sinusoids as 
shown in Figure 3.4. 
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Figure 3.4 Input Signal x(k) 
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a. ■ System I: y(k) = x(k) - (0.6y(k-l) +0.08y{k-2)) 

This is the linear system first introduced in 

Section II. D. 6. An input Gaussian white noise sequence was 

used to model the system. Figure 3.5 is a plot of the 

running average mean-square error between the input noise 

and output error signals. In Section II. D. 6, the linear 

lattice filter’s reflection coefficients were found to be 

K = -0.555 and K = -0.077. It should be noted that these 
1 2 

values appear in the first row of the nonlinear lattice’s 
reflection coefficient matrix. Here, for a first order, 
nonlinear lattice filter, the reflection coefficients are 
given by the upper triangular matrix 



K 



0 0 -.55 -.07 0 

0 0 0 0 -.43 

000 -.55 0 

0 0 0 0 0 

0 0 0 0 0 



The system output, y(k), corresponding to the input signal 
of Figure 3.4 is shown in Figure 3.6. As can be seen by the 
plots of x(k) and x(k) in Figure 3.7, the nonlinear lattice 
filter did an outstanding Job of recovering the "unknown" 
input signal x(k) from the linear system’s output y(k). 

b. System II: y(k) = x<k) - ( 0 . 2y ( k- 1 ) y ( k-2 ) ) 

This system involves a "cross-talk" 
nonlinearity, that is, the output is a function of the 
product of two different signals. For this system, both 
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Figure 3.5 System I Mean-Square Error 
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Figure 3.6 System I Output y(k) 
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Figure 3.7 System I Comparison of x(k) and x(k) 
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Gaussian and uniform noise distributions were used in the 



modeling process in order to compare the two techniques. 
Here, both methods yielded the same reflection coefficients: 



K 



0 0 
0 0 
0 0 
0 0 
0 0 



0 0 -.21 
0 0 0 

0 0 0 

0 0 0 

0 0 0 



The only difference between the two techniques was the value 
of the scaling factor (i.e., the norm of the noise generated 
system output y(k)). The plot of y(k) is shown in Figure 
3.8. The running average mean-square error and x(k) plots 
for the Gaussian noise derived model are depicted in Figures 
3.9 and 3.10, respectively. The corresponding plots for the 
uniform noise derived model are given in Figures 3.11 and 
3.12. As can be seen by comparing the plots, both modeling 
variations lead to nearly identical results. Both 
techniques yielded excellent approximations to the input 
signal X ( k ) . 

c. System III: y(k) = x(k) - 0 . 2y ( k-1 ) y ( k- 1 ) 

This system involves a quadratic nonlinearity. 
Therefore, a second order model was used. The system output 
would not remain bounded for a Gaussian noise input, so the 
system was modeled using a unit variance uniform noise 
random process. The resulting reflection coefficients are: 
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Figure 3.8 System II Output y(k) 



103 




(>|]3SW 



Figure 3.9 System II Mean-Square Error 
(Gaussian Noise Model) 
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Figure 3.10 System II Comparison of x(k) and x(k) 
(Gaussian Noise Model) 
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Figure 3.11 System II Mean-Square Error 
(Uniform Noise Model) 
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Figure 3.12 System II Comparison of x(k) and x(k) 

(Uniform Noise Model) 



107 





0 


-.20 


.08 


0 


-.04 


- . 18 


0 


0 


0 


0 




0 


0 


-.20 


-. 18 


0 


.69 


.54 


. 17 


. 34 


-.42 




0 


0 


0 


.12 


-.22 


-.39 


- .04 


.01 


.64 


.31 




0 


0 


0 


0 


-.38 


-.09 


- . 37 


. 60 


- .15 


. 32 


K = 


0 


0 


0 


0 


0 


.27 


.15 


-.67 


-.25 


-.13 




0 


0 


0 


0 


0 


0 


. 57 


-.17 


-.37 


.49 




0 


0 


0 


0 


0 


0 


0 


- . 36 


- . 37 


. 50 




0 


0 


0 


0 


0 


0 


0 


0 


.51 


-.38 




0 


0 


0 


0 


0 


0 


0 


0 


0 


-.81 




0 


0 


0 


0 


0 


0 


0 


0 


0 


0 



Figure 


3.13 


provides 


the 


running 


average 


mean-square error 


plot . 


The output y(k) 


of 


the nonlinear 


system 


is shown 


in 


Figure 


3.14. 


Figure 3 


. 15 


depicts 


the comparison 


between 


the 


system 


input 


X ( k ) and 


the 


lattice 


filter 


output 


x( k) . 


The 



x(k) curve has the same shape as x(k), however, it is 
slightly offset. In an attempt to improve the approximation 
process, the number of realizations of the noise random 
process used to model the system was doubled from twenty- 
five to fifty and the simulation was repeated. Although 
several reflection coefficient values changed by .01, there 
was no perceptible improvement in the estimation of x(k). 

d. System IV: y(k)=x(k) - ( 0 . 5+0 . 6y ( k-1 ) + . 08y ( k-2 ) ) 
This example consists of the linear System I 
modified by the addidtion of a constant bias term. Gaussian 
noise was used to model the system. The reflection 
coefficients determined for the first order model are: 
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Figure 3.13 System III Mean-Square Error 
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Figure 3.14 System III Output y(k) 
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Figure 3.15 System III Comparison of x(k) and x(k) 
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K 



0 

0 

0 

0 

0 



23 -.55 -.07 0 

0 -.23 -.40 -.43 

0 0 -.46 .01 

0 0 0 0 

0 0 0 0 



This is basically the same matrix as obtained for System I, 
but modified by the addition of several terms which attempt 
to model the constant offset. (Note the reappearance of the 
values of -.55 and -.07 in the top row of the lattice 
matrix.) The running average mean-square error is shown in 
Figure 3.16. The system output y(k) is plotted in Figure 
3.17. Figure 3.18 is a plot of the lattice output x(k) 
and the system input x(k). The x(k) curve is an excellent 
replica of the original input signal, but it is offset by a 
constant value. Although the small mean-square error evident 
in Figure 3.16 indicates that the lattice is a good inverse 
filter, the lattice was unable to completely remove the bias 
term. Increasing the order of the model (i.e., overmodeling 
the system) did not improve the results. 

e. System V: y(k)=x(k) - ( 5 . 0+0 . 6y (k-1 ) + . 08y { k-2 ) ) 

In order to further investigate the constant 
offset nonlinearity, the example of System IV was repeated 
using a bias of 5.0 instead of 0.5. Again, Guassian white 
noise was used in determining the model's parameters. The 
reflection coefficients of the first order nonlinear lattice 
are given by: 
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Figure 3.16 System IV Mean-Square Error 
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Figure 3.17 System IV Output y(k) 



114 



o 




(>niuHx'(>i]x 



Figure 3.18 System IV Comparison of x(k) and x(k) 
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K 



0 --92 -.55 -.07 0 

0 0 -.92 -.86 -.74 

0 0 0 .7-8 -.7 5 

0 0 0 0 -.90 

0 0 0 0 0 . 

Again, the reflection coefficient values of -.55 and -.07 
occur in the top row of the matrix as these elements 
apparently model the linear portion of the sytem. The 
running average mean-square error is displayed in Figure 
3.19, and the output of the nonlinear system is shown in 
Figure 3.20. Figure 3.21 provides the comparison between the 
system input and the output of the inverse filter. As with 
System IV, the x(k) curve produced by the deconvolution 
process has the same shape as x(k), but is offset from the 
desired curve by a small constant value. Here, the system 
bias is 5.0, and the x(k) and x(k) curves differ by about 
0.5. This is a relative improvement over the results 
obtained for System IV. System IV had a constant bias of 
only 0.5 but the x(k) curve was offset from the x(k) curve 
by approximately 0.4. 
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Figure 3.19 System V Mean-Square Error 
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Figure 3.20 System V Output y(k) 
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Figure 3.21 System V Comparison of x(k) and x(k) 
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IV. CONCLUSION 



In this thesis, several linear least-squares 
deconvolution, or inverse filtering, techniques were 
reviewed. Particular emphasis was placed on the lattice 
filter. It was shown that when a system is defined by an 
autoregressive, linear time-invariant model, this model can 
be transformed into an equivalent lattice filter 
representation. Furthermore, the resulting analysis lattice 
filter acts as a whitening filter, and is the inverse of the 
system^ s autoregressive transfer function. An example 
demonstrated the effectiveness of the lattice filter in a 
linear deconvolution application. Unfortunately, the linear 
lattice filter is unable to model nonlinear systems and, 
therefore, is not a viable nonlinear inverse filtering 
technique . 

The discussion of the linear lattice filter led to the 
development of the generalized lattice filter, which in turn 
led to the derivation of the nonlinear lattice filter. 
The goal of this thesis was to implement the nonlinear 
lattice filter, and then apply it to the linear and 
nonlinear deconvolution problem. This was accomplished with 
generally very good results. It was shown that the 
nonlinear lattice filter was suitable for modeling 
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discrete autoregressive systems where the output y(n) 
consisted of a weighted sum of the cross-products of 

P q 

the terms y (n-1) and y (n-2), where the powers p and q can 
take on the values {0,1, 2, 3, 4). Examples were presented 
demonstrating the ability of the nonlinear lattice filter to 
effectively act as an inverse filter for both linear and 
nonlinear autoregressive systems. 
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APPENDIX A 

LINEAR LATTICE FILTER FORTRAN PROGRAMS 






* 

* 

* 

* 



LT SCOT L. JOHNSON 
LINEAR INVERSE 
06 APRIL 1986 



*>«cyc*yc*yfyf*ycyt>«cyf>«{**ycy{>«f**ycyf***ycycycy{***yc*******ycy{yfyc****yc****y{i«fyc*yf*yf*yfycyc* 

THE PURPOSE OF THIS PROGRAM IS TO MODEL A LINEAR AUTOREGRESSIVE 
SYSTEM BY THE EQUIVALENT ANALYSIS LATTICE FILTER. 

THE system's TRANSFER FUNCTION IS GIVEN BY H(Z) = 1/A(Z) 

WHERE AfZ) = (Z**2 + A1*Z + A2)/Z**2. 

THE LATTICE FILTER REFLECTION COEFFICIENTS ARE DETERMINED 
BY DRIVING HfZ} WITH ZERO MEAN, UNITY VARIANCE GAUSSIAN WHITE 
NOISE. THE OUTPUT IS PROCESSED BY SUBROUTINE ATOCOR WHICH COMPUTES 
THE COMPONENTS OF THE SAMPLED AUTOCORRELATION MATRIX R. SUBROUTINE 
LEVIN IMPLEMENTS THE LEVINSON ALGORITHM AND EVALUATES THE LATTICE 
FILTER COEFFICIENTS FROM THE AUTOCORRELATION MATRIX. FINALLY, 
SUBROUTINE LATICE IMPLEMENTS THE ANALYSIS STRUCTURE OF THE 
LATTICE FILTER WHICH IS EQUIVALENT TO THE INVERSE FILTER H(Z). 

NOW WHEN H(Z) IS DRIVEN BY AN UNKNOWN SIGNAL XtN}. THE 
SYSTEM OUTPUT Y(N) CAN BE FED INTO THE PREVIOUSLY DETERMINED 
LATTICE FILTER: SINCE THE LATTICE FILTER IS THE INVERSE OF H(Z) 

THE LATTICE FILTER OUTPUT SHOULD BE A GOOD ESTIMATE OF X(N). 

***VARIABLE DEFINITIONS**** 

Al, A2 = COEFFICIENTS OF THE PREDICTION ERROR POLYNOMIAL A(Z) 
DS£ED,DSEED1 = SEED VALUES USED BY THE IMSL WHITE GAUSSIAN NOISE 
GENERATOR FUNCTION GGNQG 

E = VECTOR OF MEAN-SQUARED PREDICTION ERRORS E(0) ,E(1) ,. . . E(ORDER) 
ESUM = VECTOR USED IN DETERMINING AVERAGE E 
GAMMA = VECTOR OF LATTICE REFLECTION COEFFICIENTS. GAMMAfl) IS 
THE REFLECTION COEFFICIENT OF THE TTH LATTICE STAGE. 
GAMSUM = VECTOR USED IN DETERMING AVERAGE GAMMA 
GRAFP = REAL VALUE WHICH DEFINES LENGTH OF X-AXIS FOR PLOTTING 
L = LOWER-TRIANGULAR MATRIX WHOSE ROWS ARE THE REVERSE OF ALL THE 
PREDICTION ERROR FILTERS FROM ORDER ZERO TO THE HIGHEST ORDER 
LSUM = MATRIX USED IN DETERMINING THE AVERAGE L MATRIX 
MSE = MEAN-SQUARE ERROR BETWEEN X(N) AND XHAT(N) 

MSEMAX = MAXIMUM VALUE OF MSE 

MSESUM = RUNNING SUM USED IN CALCULATING MSE 

NINDEX = ARRAY OF REAL NUMBERS USED IN DISSPLA PLOTTING ROUTINES 
NOISE = ARRAY OF MEASUREMENT NOISE ADDED TO OUTPUT OF H(Z) 

NUMPTS = NUMBER OF POINTS IN INPUT NOISE SEQUENCE 
ORDER = ORDER OF THE LATTICE FILTER 
ORDERP = ORDER + 1 

PLTPTS = NUMBER OF POINTS USED IN PLOTTING ROUTINES 
R = VECTOR OF AUTOCORRELATION LAGS R(0i,R(l),. . . RfORDER) 

RMAX = MAXIMUM MAGNITUDE OF ELEMENTS OF R^ T6 BE USED IN PLOTTING 
STDDEV = STANDARD DEVIATION OF THE MEASUREMENT NOISE 
TRIAL = NUMBER OF REALIZATIONS OF THE WHITE GAUSSIAN NOISE RANDOM 
PROCESS USED IN MODELING THE SYSTEM H(Z) 

X = INPUT SEQUENCE INTO THE SYSTEM H(Z) 

XMAX, XMIN = RANGE OF X VALUES USED IN DISSPLA PLOTTING ROUTINES 
XHAT = ESTIMATE OF X: OUTPUT OF THE ANALYSIS LATTICE FILTER A(Z) 

Y = OUTPUT SEQUENCE 6F HfZ) 

YMAX, YMIN = RANGE OF Y VALUES USED IN DISSPLA PLOTTING ROUTINES 

***VARIABLE DECLARATIONS**** 

INTEGER I, N,K, NUMPTS, ORDER, ORDERP, PLTPTS, TRIAL 
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c 

c 

c 



c 

c 



4 

6 

C 



C 

c 

c 



10 



11 

c 



12 

C 

C 



18 

C 



REAL*4 XfSOOOl.YCSOOOl.Al.AZ.XMIN.XMAX.YMAXJMIN.NINDEXfSOOO) 
REAL*4 GAMMA( 16 }'gAMS 0M( l6) . 10) L( 10 10} M 10} ’ RMAX .GRAFP 
REAL*4 XHAT(5000) ,LSUM(.10 J6)'ESUrJl(10) NOI^E(5006) ,STl!)DEV 
PEAL*4 MSE(5000Kr^SESUH,M^EPLt 
REAL*8 DSE^,DSEfDl 

DEFINE THE LENGTH OF THE TIME SEQUENCES 
NUMPTS = 5000 
PLTPTS = 2000 
GRAFP = FLOAT (PLTPTS+1) 

DEFINE THE ORDER OF THE LATTICE FILTER AND THE A(Z) COEFFICIENTS 
0RDER=2 

0RDERP=0RDER+1 
A1 = 0.6 
A2 = 0.08 

INITIALIZE VARIABLES 
DO 6 I=1,0RDERP 

GAMMA(I) = 0. 

GAMSUM(I) = 0. 

ESUM(I) = 0. 

RU) - 0- 

DO 4 K=1,0RDERP 
LSUM(I,K) = 0. 

CONTINUE 

CONTINUE 

XMAX = 0. 

YMAX = 0. 

MSESUM = 0. 

MSEMAX = 0. 

DEFINE THE SEEDS FOR THE NOISE SIGNALS. THE NUMBER OF NOISE 
REALIZATIONS USED IN MODELING H(Z), ANl) THE STANDARD DEVIATION OF 
THE MEASUREMENT NOISE 
DSEED = 1243073. 5D0 
DSEEDl = 724389. 4D0 
TRIAL = 25 
STDDEV = 0. 05 

DO 25 I=1.TRIAL , ^ 

DSEEd = DSEED/DFLOAT(I) ^ 

DSEEDl = DSEED1/DFL0AT(I) 

DO 10 K=l, NUMPTS 

X(K) = GGNQF(DSEED) 

CONTINUE 

DO 11 K=L NUMPTS 

NOISd(K) = STDDEV * GGNQF(DSEEDl) 

CONTINUE 



[ 2 ] = x(2] - Al*Y(iP+ N0ISE(2) 
0 12 K = 3, NUMPTS 



Y(K) = X(K) - (A1*Y(K-1) 
CONTINUE 



A2*Y(K-2)) + NOISE(K) 



DETERMINE THE AUTOCORRELATION VECTOR OF THE SEQUENCE Y(N) 
CALL ATOCORC NUMPTS ORDER, R.RMAX) 

DETERMINE THE L MATRIX AND E VICTOR 
CALL LEVIN(ORDER,R,L,E) 

DO 18 K=l, ORDER 



CONTINUE 



GAMMMK)=-1*L(K+1,1)^ , ^ 

GAMSUM(K) = GAMSUM(K) + GAMMA(K) 



DO 22 N=1.0RDERP , ^ 

ESUM(N) = ESUM(N) + E(N) 

DO 20 K=1,0RDERP 

LSUM(N,K) = LSUM(N,K) 



L(N,K) 
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20 

22 

25 

C 

C 



26 

28 



30 

C 

C 

C 

C 

C 

C 

C 

C 

35 



40 



45 

C 

C 

48 



50 

55 

C 

C 

C 

C 



85 

C 



90 

C 

C 

C 

C 



93 

C 

C 



CONTINUE 

CONTINUE 

CONTINUE 

CALCULATE THE AVERAGE VALUES OF E, L, AND GAMMA 
DO 28 N=l,ORDERP 

E(N) = ESUM(NVFLOAT(TRIAL) 

DO 26 K=l,ORD£RP 

L(N K) = LSUM(N,K)/FLOAT(TRIAL) 

CONTINUE 
CONTINUE 
DO 30 K=l, ORDER 

GAMMA(K) = GAMSUM(K)/FLOAT(TRIAL) 

CONTINUE 

WITH INPUT Y AND LATTICE PARAMETERS GAMMA, DETERMINE THE 
OUTPUT OF THE ANALYSIS MODEL LATTICE FIlTER 
CALL LATICE(NUMPTS, ORDER, GAMMA, Y.XHAT) 

PRINT RESULTS 

PRINT THE R.E, AND GAMMA VECTORS 
WRITE(8,35) 

FORMATCtI,'!' ,T4, 'N' ,T10,'R(N)' ,T20,'E(N)' ,T30,'GAMMA(N)') 
NULL=0 

WRITE(8.40) NULL^R(l)^Ejfl) 

FORMATCtI 'O' J4 j2,Tld,3F10. 4) 

DO 45 K=2,OR0ERf5 
N=K-1 

WRITE(8,40) N,R(K),E(K),GAMMA(N) 

CONTINUE 



PRINT 
WRITE 
FORMAT, 

DO 55 1= 



THE L(I,J) MATRIX 
^{ 8 , 48 ) 

CtI 'I' 'L(I,J)=') 

j i=1,or0erp 

WRltE(8,50) {L(I.J),J=1.0RDERP) 
FORMAT(tl '0^i0CFf0.4,4Xj) 



CONTINUE 

DO THE DECONVOLUTION PROBLEM. FIRST DEFINE THE DETERMINISTIC 
INPUT X(N), THEN DETERMINE THE SYSTEM OUTPUT. 

DO 85 K=1,PLTPTS 

X(K) = 1.0 * SIN(0. 0126*FLOAT(K)) 

CONTINUE 



(1) = X(l) + NOISE(l) 

(2) = X(2) - A1*Y(1) + NOISE(2) 
0 90 K=3,RLTPTS 



Y(K) = X(K) - (A1*Y(K-1) + A2*Y(K-2)) + NOISE(K) 



CONTINUE 

INPUT Y(N1 INTO THE LATTICE TO RECOVER THE ESTIMATE OF X(N) 
CALL LAT1CE(PLTPTS, ORDER, GAMMA, Y,XHAT) 

CALCULATE THE MEAN-SQUARE ERROR (MSE) BETWEEN X(N) AND XHAT(N) 
DO 93 K=1,PLTPTS 

MSE9UM = (X(K)-XHAT(K1)*(X(K)-XHAT(K)) + MSESUM 
MSEl^KLr, SQRT(mSESUM/R 

continue'"^^ 



1SEMAX) MSEMAX = MSE(K) 



DEFINE PLOTTING PARAMETERS 
DO 95 K= l.PLTPTS 

NIND£X(K) = FLOAT(K-l) 
IF(ABS(X(K" ^ 



IF(ABS(XHA 



GT. XMAX) XMAX = ABS 

" T. xm; 



,K}).GT. XMAX) XMAX , 

IF(ABS(Y(K)).GT. YMAX) YMAX = ABS(Y(R)) 






AT(K)) 
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95 CONTINUE 

XMIN = -1.0* XMAX 
YMIN = -1.0 * YMAX 

**** DISSPLA PLOTTING ROUTINES **** 

PLOT THE SYSTEM INPUT AND THE LATTICE OUTPUT 
CALL TEK618 
CALL COMPRS 



CALL RESET('ALL') 
CALL PAGE(8.0,6. 5 
CALL XINTAX 



) 



CALL AREA2D(6.75,5. 0} 

CALL XNAMEPtIME SAMPLES$M00) 

CALL YNAME('X(N)$^,100) 

CALL CROSS ' ' 

CALL GRAECO. ,500. ,GRAFP. XMIN, ' SCALE' , XMAX) 
CALL CURVECNiNDEX X,PLT?>TS,0) 

CALL ENDPL(O) 



PLOT THE SYSTEM OUTPUT Y(N) 
CALL NOBRDR 



CALL AREA2D(6. 75,5. 0) 

CALL XNAMEPTIME SAMPLES$' ,100) 

CALL YNAMEf 'Y( n5$' ,100) 

CALL HEADINr'YCN) = X(N) - CO. 6*Y(N-l)+0. 08*YCN-2))$' .100,1. 5,1) 
CALL HEADINC ' V(N)=X( N)-(0. 6*YCN-1 )+0. 08*Y(N-2 ) )+n6is 8$ ' , lOO , 1. 5 , 
CALL HEADINC 'NOISE = N(0,O.OO25)$\lOO, 1.5, 2) 

CALL CROSS » > » / 

CALL GRAECO. ,500. ,GRAEP,YMIN, ' SCALE' ,YMAX) 

CALL CURVECNiNDEX Y,PLTf>TS,0) 

CALL ENDPL(O) 

PLOT LATTICE FILTER OUTPUT, XHAT(N) 

CALL AREA2DC6. 75,5. 0) 

CALL XNAMEPtIME SAMPLES!' ,100) 

CALL YNAME('XHAT(N)$',lO0) 

CALL GRAECO. ,500. ,GRAFP,XMIN,' SCALE' ,XMAX) 

CALL CURVECNtNDEX XHAT,f>LTPTS,0) 

CALL ENDPL(O) 

PLOT THE MEAN-SQUARE ERROR BETWEEN X(N) AND XHAT(N) 

CALL AREA2D(6.75,5. 0) 

CALL XNAMEC^TIME SAMPLES!' ,100) 

CALL YNAME('MSE(N)$' ,100) 

CALL CROSS 

CALL GRAECO. ,500. ,GRAFP,^0. , 'SCALE' .MSEMAX) 

CALL GRAECO! ,500! ,GRAFf>,o! ,' SCALE' ,0. 08) 

CALL CURVECNtNDEX MSE,PLTPtS,0) 

CALL ENDPL(O) 

CALL DONEPL 

STOP 

END 



2 ) 



SUBROUTINE ATOCOR 

GIVEN A TIME SEQUENCE YCN), THIS PROGRAM CALCULATES THE SAMPLED 
AUTOCORRELATION MATRIX TERl^S R(0) ,R( 1) ,. . . ,R(ORDER). 

WRITTEN 06 APRIL 1986 

SUBROUTINE ATOCORCNUMPTS ,Y, ORDER, R,RMAX) 

C VARIABLE DEFINITIONS: 
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0 

5 



Y(5000)= INPUT SEQUENCE 

Rho) = AUTOCORRELATION MATRIX 

NUMPfS = NUMBER OF POINTS IN THE SEQUENCE Y(N) 

ORDER = MAXIMUM LAG FOR WHICH R IS EVALUATED 

ORDERP = 0RDER+1:THE LENGTH OF THE R VECTOR 

SUM= RUNNING TOTAl OF PRODUCTS FOR A GIVEN LAG 

RMAX= MAXIMUM VALUE OF R; USED FOR DISSPLA PLOTING 

INTEGER NUMPTS, I. J,K,N, ORDER, ORDERP 
REAL Y(5000),R(lC)) SUM RMAX 
RMAX=0. 

ORDERP=ORDER+l 
DO 75 K=l, ORDERP 
SUM=0. 

N=NUMPTS-K+1 
DO 70 J=1,N 

SUM=SUM+(Y(K+J-1)*Y(J)) 

CONTINUE 



CONTINUE 

RETURN 

END 



Rj.K)=SUM/NUMPTS 



(ABS(R(K)). GT. RMAX) RMAX=ABS(R(K)) 






SUBROUTINE LEVIN 

THIS SUBROUTINE IMPLEMENTS LEVINSON'S ALGORITHM. IT GENERATES ALL 
THE PREDICTION ERROR FILTERS UP TO A GIVEN ORDER, FROM THE 
AUTOCORRELATION LAGS. 

BASED ON A PROGRAM WRITTEN BY S.J. ORFANIDIS (REF. 1: P. 333) 
MODIFIED 06 APRIL 1986 

*********************************************************** *ycyc**ycycycycy{ycyt 

SUBROUTINE LEVINf ORDER, R,L,E) 

***VARIABLE DEFINITIONS^*^* 

ORDER = ORDER OF LATTICE FILTER 

R = VECTOR OF AUTOCORRELATION LAGS Rf 0) ,Rf 1) . . . . ,R(ORDER) 

L = UNIT LOWER-TRIANGULAR MATRIX. ITS i-TH ftOW HOLDS THE I-TH 
PREDICTION-ERROR FILTER IN REVERSE ORDER. ITS FIRST COLUMN 
HOLDS THE NEGATIVES OF THE REFLECTION COEFFICIENTS, 

GAMMA(I) = -L(I + K1) FOR 1=1,2,. .. ,ORDER 
E = VECTOR OF PREDICTION ERRORS E(o) , E( 1},. . . ,E(ORDER) . 

THE MATRIX L AND THE DIAGONAL MATRIX DFOftMED BY THE E'S DEFINE 
A UL CHOLESKY FACTORIZATION OF THE INVERSE OF THE AUTOCORRELATION 
MATRIX: R INVERSE = L TRANSPOSE * D INVERSE * L 

***VARIABLE DECLARATIONS**** 

REAL R(10),E(10),L(10^10)^GAP,GAMMA 
INTEGER rlPtUS,lMINU$,J, Order, ORDERP 

orderp=orOer+i 

: SET THE UPPER TRIANGLE OF THE L MATRIX TO ZERO 

DO 60 1=1, ORDER 

iplOs=i+i 

IMINUS=I-1 

DO 60 J=IPLUS.ORDERP 
L(I,J)=0. 

CONTINUE 

L(l,l)=l. 

L(2 2)=1.0 



60 

C 
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63 



64 



68 

C 



GAP=0. 

IMINUS=I-1 
DO 63 K=1,IMINUS 

GAP=GAP+R(K+1)*L(I-1,K) 

CONTINUE 

GAMMA=GAP/E(I-1) 



L(I,1)=-1*GAMMA 
DO 64 K=2,IMINUS 



CONTINUE 



l(i;k)=l(i-i,k-i)-gamma*l(i-i,i-k) 



CONTINUE 

RETURN 

END 



L(I,n=1.0 

E(l)=E(I-l)*(l-GAMMA**2) 



*********************************************************************** 



SUBROUTINE LATICE 

THIS PROGRAM IMPLEMENTS A SINGLE CHANNEL LATTICE STRUCTURE. 

WHEN GIVEN THE LATTICE COEFFICIENTS AND THE INPUT SEQUENCE, 

THE PROGRAM DETERMINES THE OUTPUT SEQUENCE. 

WRITTEN 06 APRIL 1986 

*********************************************************************** 

SUBROUTINE LATICEfNUMPTS, ORDER, GAMMA, Y, OUTPUT) 

***VARIABLE DEFINITIONS**’^* 

NUMPTS= NUMBER OF POINTS IN THE SEQUENCES; MAX IS 5000 
ORDER= ORDER OF THE LATTICE: MAX IS 9 
GAMMA(ORDER)= LATTICE COEFFICENT ARRAY 
F = FORWARD ERROR 
B = BACKWARD ERROR 

DELAY?ORDER)= ARRAY OF DELAYED BACKWARD ERROR SIGNALS 
TEMP(ORDER) = ARRAY WHICH TEMPORARILY HOLDS THE BACKWARD ERROR 
Y(NUMPTS)=INPUT DATA ARRAY 
OOTPUT(N(JMPTS)=ARRAY OF LATTICE OUTPUT DATA 

***VARIABLE DECLARATIONS**** 

INTEGER NUMPTS,ORDER.I,K,M ^ .... 

REAL GAMMAflO) F,B,d6lAY(10),TEMP(10),Y(5000),OUTPUT(5000) 

: INITIALIZE ArRAVs 

DO 80 1=1, ORDER 
DELAY(I)=0. 

TEMP(I)=0. 

CONTINUE 

DO TIME ITERATION 
DO 88 K=^ NUMPTS 

FOR~EACH TIME INSTANT, RECURSIVELY INCREASE THE LATTICE ORDER 
DO 85 M=KORDER 

tem^M=b. 

f=f-[ga^iMAC!^I*delay(m^ ; 



80 

c 

c 



85 

88 



CONTINUE 
OUTPUT(K)=F 
CONTINUE 
RETURN 
END 



CM)-(GAMMA(M)*F) 
. _ MMAfM}*DELAY(M)) 
DELAY(M)=TEMP(M) 
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APPENDIX B 

NONLINEAR LATTICE FILTER FORTRAN PROGRAMS 



* 

PROGRAM NLMAIN * 

* 

mis PROGRAM DEFINES THE SYSTEM PARAMETERS AND THE SYSTEM'S 
[NPUTS AND OUTPUTS. IT CALLS SUBROUTINES TO DETERMINE THE * 

THE CORRESPONDING NONLINEAR ANALYSIS LATTICE MODEL. 



WRITTEN 30 APRIL 1986 



yc 
* 
* 
* 

THIS PROGRAM INPUTS A WHITE NOISE SEQUENCE X(K) INTO A AUTOREGRES- 
SIVE^ NONLINEAR SYSTEM HfZ) . THE SYSTEM^S OUTPUT, Y(K),IS 
PROCESSED TO DETERMINE THE AUTOCORRELATION MATRIX AND THE NONLINEA 
LATTICE'S REFLECTION COEFFICIENT MATRIX. SINCE THE OUTPUT OF THE 
LATTICE IS WHITE NOISE, THE LATTICE IS EQUIVALENT TO THE INVERSE 
OF THE SYSTEM. 



ONS**** 

HE IMSL WHITE GAUSSIAN NOISE FUNCTION GGNQF 



*** VARIABLE DECLARAT 
DSEED = SEED USED BY 

GRAFP = MAX X-AXIS VALUE FOR X AND Y GRAPHS 
GRAFN= MAX X-AXIS VALUE FOR MSE GRAPH 

H = ARRAY OF AUTOREGRESSIVE PARAMETERS THAT DEFINE THE SYSTEM 
lY = SEED USED BY THE UNIFORM WHITE NOISE FUNCTION URAND 
MN = N * N 

MNPl = MN + 1; NUMBER OF ROWS IN THE NONLINEAR LATTICE 
MSE = MEAN SQOARE ERROR BETWEEN THE INPUT NOISE SIGNAL AND THE 
FORWARD ERROR SIGNAL FROM THE TOP ROW OF THE LATTICE 
MSESUM = RUNNING TOTAL USED IN CALCULATION MSE 
MSEMAX = MAXIMAUM VALUE OF MSE FOR USE IN PLOTTING 
MSEPLT = REAL*4 VALUES OF MSE USED IN DISSPLA PLOTTING 
NORM = ARRAY OF FACTORS THAT NORMALIZE THE LATTICE INPUTS 
N = DIMENSION OF THE SQUARE Y TENSOR MATRIX 
NUMPTS = NUMBER OF POINTS USED IN CALCULATING THE RHO MATRICES 
NINDEX = TIME INDEX ARRAY USED IN DISSPLA PLOTS 
PLTPTS = NUMBER OF DATA POINTS USED IN DISSPLA PLOTS 
RANGE = +/- RANGE OF UNIFORM WHITE NOISE 
RHO = ARRAY OF REFLECTION COEFFICIENTS 
RHOSUM = ARRAY USED IN CALCULATING THE AVERAGE RHO MATRIX 
SCALE = RECIPROCAL OF THE NORM OF THE SYSTEM OUTPUT FOR WHEN 
THE SYSTEM IS EXCITED BY WHITE NOISE 
STDDEV = STANDARD DEVIATION OF GUASSIAN WHITE NOISE 
TRIAL = NUMBER OF NOISE REALIZATIONS USED IN CALCULATING AVG RHO 
X = INPUT INTO SYSTEM H(Z) 

XHAT = LATTICE OUTPUT THAT APPROXIMATES X 

XHPLOT = ARRAY OF REAL*4 VALUES OF XHAT USED IN DISSPLA PLOTS 
XMAX^XMIN = RANGE OF X AND XHPLOT VALUES USED IN DISSPLA PLOTS 

Y = Output of hcz} 

YPLOT = ARRAY OF REAL*4 VALUES OF Y USED IN DISSPLA PLOTS 
YMAX,YMIN = RANGE OF YPLOT VALUES USED IN DISSPLA PLOTS 
YSUM = RUNNING TOTAL USED IN CALCULATING THE SYSTEM OUTPUT 

***VARIABLE DECLARATIONS**** 

INTEGER TRIAL, NUMPTS, PLTPTS, lY.IMl.JMl 

REAL*4 X(5000), MSEMAX, RANGE ’ST0DEV!YMAX,YMIN,XMAX.XMIN,GRAFN 
REAL*4 XHPLOT($000) ,NiNDEX(^000},M^EPLT(5000} . YPLOt(50(!)0) , GRAFP 



PEAL*8 Yf 5000). R(26:26},RH0SUM(26, 26), XHAT( 5000). N0RM(26), MSESUM 
REAL*8 RHO(26,26),M^E(5000),DS£ED ALRHA(26,26),B£TA(26,26),H(5,5) 
REAL*8 SCALE, YSUM ^ 
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SET THE NUMBER OF WHITE NOISE REALIZATIONS USED TO CALCULATE 
THE AVERAGE RHO MATRIX 
TRIAL = 25 

DEFINE MODEL PARAMETERS 
N = 2 

MN = N * N 
MNPl = MN + 1 
DEFINE H(Z) COEFFICIENTS 
THE AR PAF ‘ ' 



AR PARAMETERS H(2,l)=0. 6 AND Hfl,2)=. 
GAMMA(2)=-. 55 AND GAMMA(3)=-. OB 



08 CORRESPOND TO 



DO 7 1=1, N 
DO 6 J=1,N 

H(1,J) = 0. 

CONTINUE 
CONTINUE 
H(1 J) = 0. 5 
H(^J) = 0.6 



H 



=^ 0^08 
H(2,^) = 0.2 
H(2 3) = 0. 5 
HO.l; = 0.2 

INITIALIZE MSESUM, MSEMAX AND RHOSUM MATRIX 
DO 2 1=1, MNPl 
DO 1 J=1,MNP1 

RH6SUM(I,J) = 0. 

CONTINUE 
CONTINUE 
MSESUM = 0. 

MSEMAX = 0. 

XMAX = 0. 

YMAX = 0. 

DEFINE SEED VALUES^ SATURATION LIMIT^ RANGE OF UNIFORM NOISE^ 
STD DEV OF GAUSS NOISE, AND NUMBER 08 POINTS IN TIME SEQUENCES 
lY = 1354 

DSEED = 1243073. 5D0 
SAT = 0. 7 
RANGE = 1. 73205 
STDDEV =1.0 
NUMPTS = 5000 
PLTPTS = 5000 
GRAFN = FLOATf NUMPTS 
GRAFP = FLOATf PLTPTS 






WRITE HEADER FOR UNIFORM NOISE INPUT 
WRITE(8j4) RANGE 

F0RMATCt1,'^INPUT white UNIFORM NOISE HAS ZERO MEAN AND RANGE OF 
*+/- * .F6. 3) 

WRITEjS^ 5) TRIAL, NUMPTS^ lY 

FORMATCTl/RHO IS AVERAGED OVER ',13,' TRIALS OF ',15,' POINTS. 

* INITIAt. SEED=^,I6) 

WRITE HEADER FOR GUASSIAN NOISE INPUT 
WRITEf8.4) STDDEV 

4 FORMAT(tl^ INPUT WHITE GAUSS NOISE HAS ZERO MEAN AND STD DEV OF', 
*F6. 3) 

WRITEf.8^5) TRIAL^NUMPTS.DSEED 

5 FORMATai/RHO I^ AVERAGED OVER ',13,' TRIALS OF ',15,' POINTS. 

* INITIAL SEED=^F10. 1) 

PRINT H MATRIX 

WRITE! 8.8) , , , 

8 F0RMATCt1/y(K)=X(K)-(TENS0R PRODUCT H*Y) WHERE H(I,J) = ') 

DO 10 1=1, N 

WRITE(8,9) (H(I,J),J=1,N) 
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9 

10 
C 
C 
C 
C 



C 

C 

C 

C 

11 

c 



F0RMAT(T10,5(F6.3,4X)) 

CONTINUE 



13 

14 

15 
C 
C 
C 



40 

45 

50 

C 



53 

54 
C 

55 



58 

60 

C 

C 

C 

C 



C 

C 

C 

C 



RUN TRIALS FOR DIFFERENT SEED VALUES TO GET AVERAGE RHO MATRIX 

DO 50 L=1,TRIAL 

DSEED = DSEED/DFLOAT(L) 
lY = lY/L 

SELECT EITHER A GUASSIAN OR UNIFORM INPUT NOISE SEQUENCE 
DO 11 K=1,NUMPTS 

THE INPUT NOISE IS UNIFORM ON (-RANGE, RANGE) WITH 
MEAN = ” 



VARIANCE=STDDEV**2 



EAN = 0 AND VARIANCE =( (2*RANGE)**2)/12) 

X(K) = 2. 0 * (URANDflYT - . 51 RANGE 
THE INPUT NOISE IS GAUSSIAN, MEAN=0 AND V 
X(K) = GGNQF(DSEED) * STDDEV 
CONTINUE 

Y|^2|^=jBL||2|j - + H(2,l)-Y(l) + H(3,l)n(l)n(l)) 

^ YSUM = 0. ^ 

DO 14 1=1, N 

DO l3 J=1,N 
IMl = I-l 
JMl = J-1 

YSUM=H(I,J)*C00RD(Y(K-1),IM1)*C00RD(Y(K-2),JM1)+YSUM 

CONTINUE ^ 

Y(K) = DBLE(X(K)) - YSUM 
CONTINUE 

***CALL SUBROUTINES**** 

DETERMINE AUTOCORRELATION MATRIX FOR Y SEQUENCE 
CALL NLCLATCY.NUMPTS.R.N) 

DETERMINE REFLECTION EaUTORS FROM AUTOCORRELATION MATRIX 
CALL SCHURf RHO^ R , ALPHA^ BETA . MNP 1 ) 

ADD TOGETHER THE EhO MATRICEE FROM EACH TRIAL 
DO 45 I=1,MNP1 
DO 40 J=KMNP1 

RHO$UM(I,J) = RHO(I,J) + RHOSUM(I,J) 

CONTINUE 

CONTINUE 

CONTINUE 

CALCULATE AVERAGE RHO MATRIX AND TRUNCATE TO TWO DECIMALS 
DO 54 I=1,MNP1 
DO 53 J=1,MNP1 



CONTINUE 
CONTINUE 

WRITE(8.55) 
F0RMATCT1,'RH0(I,J) 
DO 60 I = 1,MNP1 



RH0(I,J) = RHOSUM(I,J)/DFLOAT(TRIAL) 
RHO^IJ) = DINT(RH0(I,J)*100. )/100. 



= ') 



WRITEC8,58lCRH0(I,J),J=l,MNPl) 
F0RMAT(T1,10F6.2) 



CONTINUE 



*** NORMALIZE THE LATTICE INPUT SIGNALS AND PASS THE *** 
NOISE GENERATED DATA THROUGH THE LATTICE FILTER. 

CALL NORMS(Y,N,NUMPTS,NORM) 

SCALE = 1.0/NOftM(l) 

CALL NLLAT(Y,RHO>,NUMPTS,NORM,XHAT) 

***EXAMINE THE WHITENING EFFECT OF THE UTTICE FILTER BY *>' 

CALCULATING THE MEAN-SQARE ERROR BETWEEN THE INPUT WHITE 
NOISE AND THE FORWARD ERROR SIGNAL OUTPUT FROM THE TOP 
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63 

C 

C 

C 

C 

C 



85 

C 



88 

89 

90 
C 
C 
C 
C 
C 
C 
C 



98 



99 

C 

C 

C 

C 

C 



ROW OF THE LATTICE. 

DO 63 K=1,NUMPTS 

MSESUM = (DBLE(XfK))-XHAT(K))*(DBLE(X(K))-XHAT(K)) + 
MSEfK) =DSQRTrMSESDM/DFLOAT([<)) 

MSEkf(K) = SHGLCMSE(K)) 

IFrMSEkT(K).GT.HSEMAX) MSEMAX = MSEPLT(K) 

CONTINUE 

NOW THAT THE INVERSE FILTER HAS BEEN EVALUATED, DO THE 
DECONVOLUTION PROBLEM. FIRST, DEFINE THE "UNKNOWN" INPUT 
SEQUENCE X(K) AND THEN GENERATE THE SYSTEM OUTPUT Y(K). 

.,PLTPTS 

LE. 250) X(K) = FL0AT(K)/250. 

LE. 750. AND. K. GT. 250) XCK) = 2. 0 - FL0AT(K)/250. 



MSESUM 



DO 85 
IF( 
IF( 

IF 

IF 

IF 

IF( 

IF 

IF( 

IF 

IF 

CONTINUi 



=1 

K. 

K. 

K. 

K. 

K. 

K. 

K. 

K. 

K. 

K. 

K. 



LE. 1250. AND. K.GT. 750) 
LE. 1750. AND. K. GT. 1250 
LE. 2000. AND. K. GT. 1750 
LE. 2400. AND. K. GT. 2000 
LE. 2700. AND. K. GT. 2400 
LE. 2900. AND. K. GT. 2700 
LE. 3000. AND. K. GT. 2900 
LE. 3500. AND. K. GT. 3000 
LE. 5000. AND. K. GT. 3500 



ul) = FL0AT(K)/250. - 4.0 
X(K) = 1.0 - FLOAT(K-1250)/250. 
xU) = FLOAT(K-1750)/25O. - 1.0 
xU) = 0. 5 
■ XU) = -0.5 
X(K) = 0.5 
X(K) = -0. 5 






0. 7*SIN(0. 0126*FLOAT(K-3001)) 
0. 9*SIN(0. 0021*FL0ATU-3501)) 



Y^2^^=jBLp^2p - (^hll) + H(2,1)*Y(1) + H(3,1)*Y(1)*Y(1)) 

^ YSUM = 0. 

DO 89 I=1,N 

DO ^8 J=1,N 
IMl = I-l 
JMl = J-1 

YSUM=H(I ,J)*C00RD(Y(K-1) ,IMl)*C00RD(Y(K-2) ,JM1)+YSUM 
CONTINUE 
CONTINUE 

Y(K) = DBLE(X(K)) - YSUM 
CONTINUE 

NOW PASS THE SYSTEM OUTPUT DATA THROUGH THE LATTICE FILTER 
EMBEDDED WITH THE PREVIOUSLY CALCULATED REFLECTION FACTORS ^ 
RHO(I^J). THE FORWARD ERROR SIGNAL OUT OF THE TOP ROW OF THE 
LATTICE IS XHATCK). AND SHOULD APPROXIMATE THE DESIRED SIGNAL 
X(K) AFTER IT IS RENORMALIZED AND SCALED. 

CALL NORMS(Y,N,PLTPTS,NORM) 

CALL NLLATf Y>HO , N , P Lt PTS , NORM , XHAT ) 

DO 98 K=l,ktPTS ^ ^ , ,, 

XHAT(I^) = SCALE * CXHATf K)*NORM( 1)) 



XHPLOT(K) = SNGL(XHAT(K 
YPLOT(K) = SNGL(Y(K)) 

IF(ABS(X(K)).GT.XMAX) XI 

V(ABS(XHPLOT(K)).GT.XM,..., 

FCABS(YPLOT(K)).GT.YMAX) YMAX = ABS(YPL0T(K 



, XMAX = ABS(X(K 

i).GT.XMAX) XMAX = AB 



XHPLOT 



C0NTINU_ 

XMIN = -1.0 * XMAX 
YMIN = -1.0* YMAX 
DO 99 K=l,5000 

NINDE)((K) = FLOAT(K-l) 
CONTINUE 






* 



DISSPLA PLOTTING ROUTINES 



PLOT SYSTEM INPUT AND LATTICE OUTPUT 
CALL TEK618 
CALL COMPRS 
CALL RESETCALL'^ 



CALL PAGE(8. 0,6. 
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c 

c 



CALL X INTAX 

CALL AREA2D(6. 75,5. 0) 

CALL XNAMEf'TIME SAMPLES!' ,100) 

CALL YNAME('X(K),XHAT(K)$' 100) 

CALL YNAMEt^X(K)$^100) 

CALL CROSS 

CALL GRAFfo. ,500, ,GRAFP ,XMIN ,' SCALE' ,XMAX) 

CALL LINESP (2. OV 

CALL LINES! 'X(K)!',IPAK.l) 

CALL LINES! 'XflAT(K)$' ,IPAK,2) 

CALL LEGLIN 
CALL DOT 

CALL CURVE!NINDEX,X,PLTPTS,0) 

CALL RESET! 'DOT') 

CALL CURVE!NINDEX,XHPLOT.PLTPTS,0) 

CALL LEGEND (IRAK 2,5. 0 , 0 . 5) 

CALL ENDPL(O) 

PLOT Y 

CALL AREA2D! 6. 75,5.0) 

CALL XNAMEPtIME SAMPLES!' ,100) 

CALL YNAME!'Y!K)$^, 100) 



C 

C 



CALL HEADINC'YrK) = X(K) - 0. 2*Y(K-1)**2$' ,100,1. 5.1) 
CALL HEADIN!'V(K) = X(K>!0. 6*Y!K-l)+0. 08’‘Y!K-2))V^ 
CALL HEADIN('^Y(K)=X(K)-(0. i*Y(K-i)+ 0. 5*Y(K-i)*Y(K-2)**2 






C 

c 



*100,1. 5,1) 

CALL CROSS 

CALL GRAF!0. ,500. ,GRAFP ,YMIN.' SCALE' ,YMAX) 

CALL CURVE!NInDEX YPLOT PLTPtS,0) 

CALL ENDPL(O) 

PLOT MSE 

CALL AREA2D(6.75,5. 0) 

CALL XNAMEPtIME SAMPLES!' ,100) 

CALL YNAME!' MSE! K)!' ,100) 

CALL GRAF!0. ,500. ,GRAFN.0. ,' SCALE' ,MSEMAX) 

CALL CURVE!n 1NDEX MSEPLt,NUMPTS,0) 

CALL ENDPL!0) 

CALL DONEPL 

STOP 

END 

SUBROUTINE SCHUR 

CALCULATES THE REFLECTION FACTORS FROM THE CORRELATION MATRIX 
WRITTEN 29 APRIL 1985 BY P.J. LENK (REF 12: P. 174) 

SUBROUTINE SCHUR(RHO,R, ALPHA, BETA.N) 

REAL*8 RH0(26,26),R(^6;26),Af.PHA(^6,26),BETA(26,26),RN0RM,T 

INITIALIZE THE ALPHA AND BETA ARRAYS 

DO 10 I = 1,N 
DO 5 J = LN 



7 

0 CONTINUE 



RHO(l,J) = 0.0 
CONTINUE 



BEGIN CALCULATING THE REFLECTION FACTORS 
DO 50 J = 2,N 
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NJl = N - J + 1 
DO 40 I = 1,NJ1 
JIl = J + I - 1 
IPl = I + 1 

RHOCI.JIl) = ALPHA(I,JI1)/BETA(IP1,JI1) 
RNORM = DSQRT(1.0 - I^H0(I,JIl)^RH0{l,JIl)) 
DO 30 K = I.N 



30 

40 

C 

C42 

C 

50 

C 

C 



CONTINUE 
CONTINUE 



T = ALP^iA(I,K) 

ALPHA(I,K) = (ALPHA('I,K)-RHOfI,JIl)*BETA(IPl,K))/RNORM 
BETA(l,ll) = (BETA(lk K)-RH0(I,JI1)*T)/RN0RM 



CONTINUE 



WRITE(8,42)J.((ALPHA(I,K),K-1,N),I=1, 

FORMATf;2X;i5,4(:2X,E12.5)) 

WRITE(8,42)J,((0ETA(I,K),K=1,N),I=1,N 



=1,N) 
N) 



RETURN 
END 

FUNCTION URAND 

TAKEN FROM "COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS" BY 
G. E. FORSYTHE, M. A. MALCOLM, AND C. B. MOLER 

************ ******>'C***)lf******2«C*****>'C*5'C5«C**>IC5«C*)'C*5'C5«C5<C*>'f)‘Olc*>IC*5«C>'C*>'C*****5«f>'C5*C*)'C* 

REAL FUNCTION URAND(IY) 

INTEGER IA,IC,ITW0,H2,M,MIC 
DOUBLE PRECISION HALFM 
REAL S 

DOUBLE PRECISION DATAN,DSQRT 
DATA M2/0/.ITW0/2/ 

IF(M2. NE.O) GO TO 20 

IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH 

M=1 
10 M2=M 

M=ITW0*M2 

IF(M.GT.M2) GO TO 10 

haLfm=m2 

COMPUT MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD 

IA=8*IDINT(HALFM*DATAN(1. D0)/8. D0)+5 
IC=2*IDINT(HALFM*(. 5D0-DSQRT(3. D0)/6. D0))+1 
MIC=(M2-IC)+M2 

S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT 
S=. 5/HALFM 

COMPUTE NEXT RANDOM NUMBER 
20 IY=IY*IA 

THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW 
INTEGER OVERFLOW ON ADDITION 

IF(IY.GT.MIC) IY=(IY-M2)-M2 

IY=IY+IC 



THE FOLLOWING IS FOR COMPUTERS FOR WHICH THE WORD LENGTH 
FOR ADDITION IS GREATER THAN FOR MULTIPLICATION 
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IF(IY/2. GT. M2)IY=(IY-M2)-M2 

C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER OVERFLOW 
C AFFECTS THE SIGN BIT 
C 

IF(IY. LT.0)IY=(IY+M2)+M2 

URAND=FL0AT(IY)*S 

RETURN 

END 

C * 

C SUBROUTINE NLCLAT 

C 

C THIS SUBROUTINE PRODUCES A CORRELATION MATRIX FROM NONLINEAR * 

C TIME SEQUENCE IN AN ORDER WHICH IS COMPATIBLE WITH SUBROUTINE * 

C SCHUR. * 

C * 

C WRITTEN 7 MAY 1985 BY P.J. LENK (REF 12: P. 181) * 

SUBROUTINE NLCLAT (YJYS,R,N) 

REAL*8 Y(5000),R(26,^6),3UM,VEC(26) 



C 

C 

C 



10 

20 

C 

C 

C 



30 

40 

50 

C 

C 

C 



C 

C12 

60 

70 

80 

C 



DEFINE CONSTANTS 

MN = N*N 
MNPl = MN + 1 
IYSM2 = lYS - 2 
FIYSM2 = FL0AT(IYSM2) 

INITIALIZE R MATRIX TO ZERO 

DO 20 I = 1,MNP1 
DO 10 J = 1,MNP1 



CONTIN 

CONTINUE 






= 0.0 



BEGIN OUTER LOOP 

DO 80 I = 3,IYS 
IR = 1 

VECfIR) = Y(I) 

DO 50 MPl = 1,N 
MO = MPl - 1 
LLIM = 2*MP1 - 1 
DO 40 L = l.LLIM 
LO = L - 1 
II = MO 
J1 = LO/2 

IF (M0D(L0,2).EQ. 0) GO TO 30 

J1 = MO 
IR = IR + 1 

VEC(IR) = C00RD(Y(I-1),I1)*C00RD(Y(I-2),J1) 
CONTINUE 
CONTINUE 

CALCULATE THE CORRELATIONS 

DO 70 J = l.MNPl 
DO 60 K = J.MNPl 

R(J,K) = RCJ.K) + VEC(J)*VEC(K) 
WRITE(6'12)VEC(J},VEC(K),R(J,K) 
F0RMAT(5(2X,E12. 5)) 

CONTINUE 

CONTINUE 

CONTINUE 
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90 

100 

C 

C 

C 



no 

120 

c 

c 



DIVIDE BY THE NUMBER OF DATA ELEMENTS CONSIDERED 

DO 100 J = l.MNPl 
DO 90 K = J.MNPl 

R(J,K) = R(J,K)/FIYSM2 
CONTINUE 
CONTINUE 

FILL IN THE SYMMETRIC HALF OF CORRELATION MATRIX 

DO 120 I = 2,MNP1 
IMl =1-1 
DO no J = 1,IM1 
R(I,J) = R(J,I) 

CONTINUE 

CONTINUE 



RETURN 

END 



C 

c 
c 

C FUNCTION COORD 

C 

C GENERATES OUTPUT OF RANDOM FUNCTION 

C CREATED 23 AUG 84 (REF 12: P. 187) 

C 

DOUBLE PRECISION FUNCTION COORD(XjI) 

USE SIMPLE POWER SERIES TYPE POLYNdMIALS 



C 

C 



30 



Y = 1.0 

IF (LEO. 0) GO TO 30 

Y = X**I 
COORD = DBLE(Y) 

RETURN 

END 



C 
C 

c 

C SUBROUTINE NLLAT 

C 

C THIS SUBROUTINE IMPLEMENTS THE NONLINEAR LATTICE FILTER USING 

C PREVIOUSLY CALCULATED REFLECTION FACTORS. 

C 

C WRITTEN 30 APRIL 86 

C 



C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



SUBROUTINE N LLATf Y . RHO, N , NUMPTS , NORM , XHAT) 

***VARIABLE DEFINITfONS’^*^* 

INFWD(I) = FORWARD ERROR INPUT INTO THE I'TH ROW OF THE LATTICE 
INBKDfl) = BACKWARD *' '■ " '* " '' '' '' " 

OUTFWD(I)= FORWARD ERROR OUTPUT FROM THE I'TH ROW OF THE LATTICE 
0UTBKDfl1= BACKWARD " " " '' " " '' 

Y = INPUT DATA VECTOR 
RHO = REFLECTION FACTOR MATRIX 
NORM = VECTOR OF NORMS OF THE LATTICE INPUT TERMS 
XHAT = OUTPUT DATA VECTOR; IT IS THE FORWARD ERROR SIGNAL FROM 
THE LAST STAGE OF THE FIRST ROW 
NUMPTS = NUMBER OF POINTS IN THE INPUT/OUTPUT SEQUENCES 
N = DIMENSION OF SQUARE Y DATA MATRIX 

MNPl = DIMENSION OF THE RHO, NORM, AND INPUT/OUTPUT ARRAYS 

***VARIABLE DECLARATIONS**** 

INTEGER N,MN, NUMPTS, LAST, MNPl 
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c 

c 



5 

C 

C 



10 



11 

12 

13 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 



14 
C 

c 

15 



C 

C 

20 



REAL*8 Y(5000)^RHO(26,26)^NORMf26},XHAT(5000),INFWD(26) 
REAL*8 INBKD(2S),OOTFWD(2^),OUTBKD(26),RNORM 

MN = N * N 
MNPl = MN + 1 

INITIALIZE THE INPUT AND OUTPUT ARRAYS TO ZERO 
DO 5 1=1.26 

INFWCf(I) = 0. 

INBKDa) = 0. 

OUTFWD(I) = 0. 

outbkdU) = 0 . 

CONTINUE 



DETERMINE THE LATTICE INPUTS FOR EACH TIME K. 
DO 60 K=1,NUMPTS 

CLEAR If^PUTS FOR THIS TIME ITERATION 
DO 10 1=1, MNPl 



CONTINUE 
IF(K. E( 
IF(K. E( 
SET LATT 
IR = 1 



INFWD(I) = 

INBKDa) = 

I. 1 ) GO TO 
1.2) GO TO 
;CE INPUTS 



0 . 

0 . 

15 

20 

FOR CASE WHEN 



K>=3 



INFWD(IRJ = Y(K)/NORM(IR) 



DO 



13 MPl = IN 
MO = MPl - 1 
LLIM = 2*MP1 - 1 
DO 12 L = l.LLIM 
LO = L - 1 
II = MO 
J1 = LO/2 

IF (MOD(LO,2). EQ. 0) GO TO 11 
II = J1 
J1 = MO 

IR = IR + 1 , , . 

INFWD(IR)=C00RD(Y(K-l),Il)*C00RD(Y(K-2) 

CONTINUE 



CONTI 
‘NFWD 
NFWD 
NFWD 
NFWD 
NFV/D 
NFWD 
NFWD 
NFWD 
NFV/D 
NFWD 
DO 

INBKD(J) 

CONTINUE 
GO TO 25 



= Y(K)/N0RM(1) 
= i:/N0RM(2) 




/NORM 
/NORM 
/NORM,^,, 

*Y(K-2)/N0RM(81 
*Y(K-2l/N0RM(9) 
)*Y(K-2)*Y(K-2)/NORM(10) 



= INFWD(J+1) 



SET INPUTS FOR K=1 CASE 
INFWD(l) = Y(1)/N0RM(1) 
INFWDC2) = l./NORMa) 
INBKDai = INFWD(2) 

GO TO 25 



SET 
INFWDfl 
I NFWD 
I NFWD 
I NFWD 
I NFWD 
I NFWD 
INBKD 
INBKD 



INPUTS F^R.K=2_CASE 



2)/N0RM(l) 

2) = l./N0RMr2) 

3) = Y(l)/N0RMf3 

1 7 Z 




1) = INFWD(2 

2) = INFWDa 



.(n/N^RM 
)*Y(1)*Y(1 

<*Y(l)*Ya 



/^ORMfin , , 

*Y(1)/N0RM(18) 



J1)/N0RM(IR) 



136 



c 

c 

c 

25 



C 

C 



35 

C 

40 



45 

50 

C. 

60 



(5), 


= INFWD(6) 


flO 


) = INFWD(ll) 


17 


) = INFWDUS) 



INBKD 



IMPLEMENT THE J'TH UPPER DIAGONAL OF THE LATTICE. MOVE ALONG THE 
DIAGONAL BY STARTING AT ROW L=1 AND MOVING DOWN TO ROW L=LAST. 

DO 50 0=1, MN 

iJST = MNPl - 0 
IFfO. EO. 1) GO TO 40 

IF NOT ON THE FIRST UPPER DIAGONAL (I.E. J NOT = 1) THEN UPDATE 



THE INPUTS FROM THE OUTPUTS OF THE PREVIOUS DIAGONAL 
DO 35 L=1,UST 

INFWD(L) = OUTFWD(L) 

INBKOa; = OUTBKOa+1) 

CONTINUE ^ 

DO THE LATTICE CALCULATIONS 
DO 45 L=1.LAST 
f^NORM = 



OSQRn I.^Od” RH0( K wrnORM 

(L) = aNBKD(L) - RHOa,J+L)*INFV/Da))/RNORM 



OUTFWD 
OUTBKD 

CONTINUE 
CONTINUE 

REMOVE NORMALIZATION FACTOR FROM THE OUTPUT OF THE TOP ROW 
XHAT(K) = OUTFWD(l) 

CONTINUE 

RETURN 

END 



C 
C 
C 
C 
C 

Qyc>if>'cyfycyc>'c^«c>»cyc><cyc><cyc><cyc><c>'cyt><cyc***yc><c>'fycyc**yc**yc**>ic*>'cyc>ic*>'c5*cyc><c*>>c?«c*><c*>'cycyc*yc*****yc**ycyc*yc 

c 

C SUBROUTINE NORMS 
C 

C THIS SUBROUTINE CALCULATES THE NORMS OF THE INPUTS 

C TO THE NONLINEAR LATTICE. 

C 

C WRITTEN 30 APRIL 86 

C 

SUBROUTINE NORMS(Y,N,NUMPTS,NORM) 

C 

c 

INTEGER NUMPTS,N,MN,MNP1,I,K.IR,LLIM,L0,M0,MP1,I1,J1,J,L 
REAL*8 VEC(26),N(!)RM{26),t(S060) 

MN = N * N 
MNPl = MN + 1 

C INITIALIZE VECTORS TO ZERO 

DO 10 K = l.MNPl 



VEC 



10 

C 

C 

C 



C 

C 

C 



NOR^ll(^) =*^0. 

CONTINUE 

SINCE TIME INDEX STARTS AT 1=3, INITIALIZE NORMS OF INPUTS WHICH 
ARE POWERS OF Y(I-l) TO ACCOUNt FOR TIMES 1=1 AND 1=2. 

NORM(l) = Y(l) * Y(l) + Y(2) * Y(2) 

NORmU) = 1.0 +1.0 
NORM(3) = Y(l) * Y(l) ^ ^ 

NORM^) = NORM(3) ^ NORM(3) 

NORM(ll) = NORM(6) * NORM(3) 

N0RM(18) = NORM(6) * NORM(6) 

FOR TIME L FORM A Y DATA VECTOR WHOSE COMPONENTS MATCH THE 

INPUTS TO The nonlinear lattice. 

DO 80 I = 3.NUMPTS 
IR = 1 

VEC(IR) = Y(I) 
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30 

40 

50 

C 

C 



60 

80 

C 

83 



86 

85 



DO 50 MPl = l.N 

MO = MPl - 1 
LLIM = 2*MP1 - 1 
DO 40 L = l.LLIM 
LO = - 1 



II = MO 
J1 = LO/2 

IF (MOD(LO,2). EQ. 0) GO TO 30 
II = J1 
J1 = MO 
IR = IR + 1 

VEC(IR) = C00RD(Y(I-1),I1)*C00RD(Y(I-2),J1) 



CONTINUE 
CONTINUE 

CALCULATE THE NORMS 
DO 60 K=1.MNP1 

IF(K. EO. 2) GO TO 60 
NORM(K) = NORM(K) 
CONTINUE 
CONTINUE 

WRITE(8.83) 

FORMAT(Tj 'K' JIO , ' NORM( K) ' ) 
DO 85 K=l,MNf>l 



VEC(K)*VEC(K) 



NORMCKf) = DSQRT(NORM(K)/DFLOAT(NUMPTS)) 
WRITEC8,86) K.NORMCK) 
F0RMAT(T1,I3,T8,E12. 5) 



CONTINUE 

RETURN 

END 
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